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Abstract 



o 
o 

Q I We review recent efforts to re-formulate the Einstein equations for fully relativistic numerical 

O^' simulations. The so-called numerical relativity (computational simulations in general relativity) is 

5^ I a promising research field matching with ongoing astrophysical observations such as gravitational 

W)' wave astronomy. Many trials for longterm stable and accurate simulations of binary compact 

objects have revealed that mathematically equivalent sets of evolution equations show different 
numerical stability in free evolution schemes. In this article, we first review the efforts of the com- 
. , munity, categorizing them into the following three directions: (1) modifications of the standard 

■ Arnowitt-Deser-Misner equations initiated by the Kyoto group, (2) rewriting of the evolution equa- 

tions in hyperbolic form, and (3) construction of an "asymptotically constrained" system. We next 
introduce our idea for explaining these evolution behaviors in a unified way using eigenvalue analy- 
sis of the constraint propagation equations. The modifications of (or adjustments to) the evolution 
equations change the character of constraint propagation, and several particular adjustments using 
constraints are expected to diminish the constraint-violating modes. We propose several new ad- 
justed evolution equations, and include some numerical demonstrations. We conclude by discussing 
some directions for future research. 



This article is for a part of the book Progress in Astronomy and Astrophysics (Nova Science Publ., 
2003?). Also available as gr-qc/0209111. 
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1 Overview 



1.1 Numerical Relativity 

The theory of general relativity describes the nature of the strong gravitational field. The Einstein 
equation predicts quite unexpected phenomena such as gravitational collapse, gravitational waves, 
the expanding universe and so on, which are all attractive not only for researchers but also for the 
public. The Einstein equation consists of 10 partial differential equations (elliptic and hyperbolic) 
for 10 metric components, and it is not easy to solve them for any particular situation. Over the 
decades, people have tried to study the general-relativistic world by finding its exact solutions, by 
developing approximation methods, or by simplifying the situations. Among these approaches, direct 
numerical integration of the Einstein equations can be said to be the most robust way to study the 
strong gravitational field. This research field is often called "numerical relativity" . 



Numerical Relativity 


Box 1.1 


= Necessary for unveiling the nature of strong gravity. For example: 




• gravitational waves from colliding black holes, neutron stars, supernovae, ... 




• relativistic phenomena like cosmology, active galactic nuclei, ... 




• mathematical feedback to singularity, exact solutions, chaotic behavior, ... 




• laboratory for gravitational theories, higher-dimensional models, ... 





Numerical relativity is now an essential field in gravity research. The current mainstream in 
numerical relativity is to analyze the final phase of compact binary objects (black holes and/or neutron 
stars) related to gravitational wave observations (see e.g. the conference proceedings [66]). Over the 
past decades, many groups have developed their numerical simulations by trial and error. Simulations 
require large-scale computational facilities, and long-time stable and accurate calculations. So far, we 
have achieved certain successes in simulating the coalescence of binary neutron stars (see e.g. [75]) 
and binary black holes (see e.g. [6]). However, people have still been faced with unreasonable numerical 
blow-ups at the end of simulations. 

Difficulties in accurate/stable long-term evolution were supposed to be overcome by choosing 
proper gauge conditions and boundary conditions. However, recent several numerical experiments 
show that the (standard) Arnowitt-Deser-Misner (ADM) approach [12, 80, 98] is not the best formu- 
lation for numerics, and finding a better formulation has become one of the main research topics. A 
majority of workers in the field now believe in the existence of constraint-violating modes in most of 
the formulations. Thus, the stability problem is now shedding light on the mathematical structure of 
the Einstein equations. 

The purpose of this article is to review the formulation problem in numerical relativity. Gener- 
ally speaking, there are many open issues in numerical relativity, both theoretical (mathematical or 
physical) and numerical. We list major topics in Box 1.2. More general and recent introductions to 
numerical relativity are available, e.g. by dTnverno (1996) [50], Seidel (1996/98/99) [73], Briigmann 
(2000) [25], Lehner (2001) [56], van Putten (2001) [88], and Baumgarte-Shapiro (2002) [16]. 
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Numerical Relativity open issues 


Box 1.2 


0. How to select the foliation method of space-time 

Cauchy (3 + 1), characteristic (2 + 2), or combined? 




if the foliation is (3 + 1), then • • • 




1. How to prepare the initial data 




Theoretical: Proper formulation for solving constraints? 
How to prepare realistic initial data? 
Effects of background gravitational waves? 
Connection to the post-Newtonian approximation? 

Numerical: Techniques for solving coupled elliptic equations? 
Appropriate boundary conditions? 




2. How to evolve the data 




Theoretical: Free evolution or constrained cvohition? 

Proper lormuiation tor the evolution equations; < 
Suitable slicing conditions (gauge conditions)? 

Numerical: Techniques for solving the evolution equations? 
Appropriate boundary treatments? 
Singularity excision techniques? 

Matter and shock surface treatments? 
Parallclization of the code? 


=<^<^ this review 


o. xlow LO exLiacL iiie piiybicai iiiioiiiiaLioii 




Theoretical: Gravitational wave extraction? 

Connection to other approximations? 

Numerical: Identification of black hole horizons? 
Visualization of simulations? 





1.2 Formulation Problem in Numerical Relativity: Overview 

There are several different approaches to simulating the Einstein equations. Among them the most 
robust way is to apply 3+1 (space + time) decomposition of space-time, as was first formulated by 
Arnowitt, Descr and Misner (ADM) [12] (wc call this the "original ADM" system). ^ 

If we divide the space-time into 3+1 dimensions, the Einstein equations form a constrained system: 
constraint equations and evolution equations. The system is quite similar to that of the Maxwell 
equations (Box 1.3), 

^One alternative method of space-time foliation is the so-called characteristic approach (2 + 2 space-time decomposi- 
tion). Sec reviews e.g. by d'Invorno (1996) [50], Winicour [89], Lehner (2001) [56]. Even in the 3-1-1 ADM approach, we 
concentrate the standard finite differential scheme to express numerical expression of space-time. See e.g. Brewin [23] 
for a recent progress in a lattice method. 
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The Maxwell equations : 

The evolution equations: {dt = d/dt) 



Box 1.3 



5tE 



rot B - 47rj, 



and 



9tB = -rot E 



(1.1) 



Constraint equations: 



div E = 47rp, 



and 



div B = 



(1.2) 



where people solve constraint equations on the initial data, and use evolution equations to follow the 
dynamical behaviors. 

In numerical relativity, this free-evolution approach is also the standard. This is because solving 
the constraints (non-linear elliptic equations) is numerically expensive, and because free evolution 
allows us to monitor the accuracy of numerical evolution. In black-hole treatments, recent "excision" 
techniques do not require one to impose explicit boundary conditions on the horizon, which is also 
a reason to apply free evolution scheme. As we will show in the next section, the standard ADM 
approach has two constraint equations; the Hamiltonian (or energy) and momentum constraints. 

Up to a couple of years ago, the "standard ADM" decomposition [80, 98] of the Einstein equation 
was taken as the standard formulation for numerical relativists. However, numerical simulations were 
often interrupted by unexplained blow-ups (Figure. 1). This was thought due to the lack of resolution, 
or inappropriate gauge choice, or the particular numerical scheme which was applied. However, after 
the accumulation of much experience, people have noticed the importance of the formulation of the 
evolution equations, since there are apparent differences in numerical stability although the equations 
are mathematically equivalent ^. 

At this moment, there are three major ways to obtain longer time evolutions. Of course, the ideas, 
procedures, and problems are mingled with each other. The purpose of this article is to review all 
three approaches and to introduce our idea to view them in a unified way. Table 1 is a list of references. 

(1) The first possibility is to use a modification of the ADM system developed by the Kyoto group 
[63, 64] (often cited as Shibata and Nakamura [74]) and later re-introduced by Baumgarte and 
Shapiro [15]. This is a combination of new variables, conformal decomposition, rescaling of 
the conformal factor, and replacement of terms in the evolution equation using momentum 
constraints (see §2.1). 

(2) The second direction is to rc-formulatc the Einstein equations in a first-order hyperbolic form. 
This is motivated from the expectation that the symmetric hyperbolic system has well-posed 
properties in its Cauchy treatment in many systems and also that the boundary treatment can be 
improved if we know the characteristic speed of the system. In constructing hyperbolic systems, 
the essential procedures are to adjust equations using constraints and to introduce new variables, 
normally the spatially derivativcd metric (sec §2.2). 

^The word stability is used quite different ways in the community. 

• We mean by numerical stability a numerical simulation which continues without any blow-ups and in which data 
remains on the constrained surface. 

• Mathematical stability is defined in terms of the well-posedncss in the theory of partial differential equations, such 
that the norm of the variables is bounded by the initial data. Sec oq. (2.34) and around. 

• For numerical treatments, there is also another notion of stability, the stability of finite differencing schemes. 

This means that numerical errors (truncation, round-off, etc) arc not growing by evolution, and the evaluation 
is obtained by von Neumann's analysis. Lax's equivalence theorem says that if a numerical scheme is consistent 
(converging to the original equations in its continuum limit) and stable (no error growing), then the simulation 
represents the right (converging) solution. See [30] for the Einstein equations. 
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time 




Constrained/ Surface 
(satisfies /Einstein's constraints) 



Figure 1: Origin of the problem for numerical relativists: Numerical evolutions depart from the 
constraint surface. 

(3) The third is to construct a system which is robust against the violation of constraints, such that 
the constraint surface is an attractor. The idea was first proposed as a "A-system" by Brodbeck 
et al [24] in which they introduce artificial flow to the constraint surface using a new variable 
based on the symmetric hyperbolic system (see §2.3). 

The third idea has been generalized by us as an asymptotically constrained system. The main pro- 
cedure is to adjust the evolution equations using the constraint equations [94, 95, 78]. The method 
is also applied to explain why the above approach (1) works, and also to propose alternative systems 
based on the ADM [95, 78] and BSSN [96] equations. Section 3 is devoted to explain this idea with 
an analytical tool of the eigenvalue analysis of the constraint propagation. 

We follow the notations of that of MTW[61], i.e. the signature of the space-time is ( — h -|--|-), and 
the Riemann curvature is defined as 



We use 1/ = 0, • • • , 3 and i,j = 1, • • • , 3 as space-time indices. The unit c = 1 is applied. The 
discussion is mostly to the vacuum space-time, but the inclusion of matter is straightforward. 




'-'p'- ua ' ^ aa^ i 



(1.3) 
(1.4) 
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formulations 


numerical applications 


(0) The standard ADM formulation 


ADM 


1962 Arnowitt-Deser-Misner [12, 80] 


=^ many 


(1) The BSSN formulation 


BSSN 


1987 Nakamura et al [63, 64, 74] 


1987 Nakamura et al [63, 64] 






1995 Shibata-Nakamura [74] 






^ 2002 Shibata-Uryu [75] etc 




1999 Baumgarte-Shapiro [15] 


=^ 1999 Baumgarte-Shapiro [15] 






=> 2000 Alcubierre et al [5, 7] 






=^ 2001 Alcubierre et al [6] etc 




1999 Alcubierre et al [8] 






1999 FrittcUi-Reula [42] 






2002 Laguna-Shoemaker [55] 


^ 2002 Laguna-Shoemaker [55] 


(2) The hyperbolic formulations 


BM 


1989 Bona-Masso [18, 19, 20] 


1995 Bona et al [20, 21, 22] 






=^ 1997 Alcubierre, Masso [2, 4] 




1997 Bona et al [21] 


2002 Bardeen-Buchman [17] 




1999 Arbona et al [11] 




CB-Y 


1995 Choquet-Bruhat and York [32] 


1997 Scheel et al [70] 




1995 Abrahams et al [1] 


1998 Scheel et al [71] 




1999 Anderson- York [10] 


2002 Bardeen-Buchman [17] 


FR 


1996 Frittelli-Reula [41] 


^ 2000 Hern [44] 




1996 Stewart [81] 




KST 


2001 Kidder-Scheel-Teukolsky [52] 


^ 2001 Kidder-Schccl-Teukolsky [52] 






2002 Calabrese et al [27] 






2002 Lindblom-Scheel [58] 




2002 Sarbaeh-Tigho [69] 




CFE 


1981 Friedrich[36] 


^ 1998 Frauendiener [35] 






1999 Hiibner [46] 


tetrad 


1995 vanPutten-Eardley[86] 


1997 vanPutten [87] 


Ashtekar 


1986 Ashtekar [13] 


2000 Shinkai-Yoneda [77] 




1997 Iriondo et al [48] 






1999 Yoneda-Shinkai [92, 93] 


2000 Shinkai-Yoneda [77, 94] 


(3) Asymptotically constrained formulations 


A-system to FR 


1999 Brodbeck et al [24] 


2001 Siebel-Hiibner [79] 


to Ashtekar 


1999 Shinkai-Yoneda [76] 


2001 Yoneda-Shinkai [94] 


adjusted to ADM 


1987 Detweiler [33] 


=> 2001 Yoneda-Shinkai [95] 


to ADM 


2001 Shinkai-Yoneda [95, 78] 


2002 Mexico NR Workshop [59] 


to BSSN 


2002 Yoneda-Shinkai [96] 


2002 Mexico NR Workshop [59] 






^ 2002 Yo-Baumgarte-Shapiro [90] 



Table 1: References to recent efforts of reformulating the Einstein equations. We list mainly those 
that have been applied to actual numerical comparisons. 
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2 The standard way and the three other roads 



2.0 Strategy 0: The ADM formulation 
2.0.1 The original ADM formulation 

The idea of space-time evolution was first formulated by Arnowitt, Deser, and Misner (ADM) [12]. The 
formulation was first motivated by a desire to construct a canonical framework in general relativity, 
but it also gave the community to the fundamental idea of time evolution of space and time: such as 
foliations of 3-dimensional hypersurface (Figure 2). This original ADM formulation was translated to 
numerical relativists by Smarr [80] and York [98] in late 70s, with slightly different notations. We refer 
to the latter as the standard ADM formulation since this version is the starting point of the discussion. 
The story begins by decomposing 4-dimensional space-time into 3 plus 1. The metric is expressed 

by 

ds'^ = g^udxV = -a^dt^ + -fijidx' + (3'dt){dx^ + (3^dt), (2.1) 

where a and Pj are defined as a = 1/ \/—g^ and Pj = goj , and called the lapse function and shift 
vector, respectively. The projection operator or the intrinsic 3-metric gij is defined as 7^;^ = g^u+n^riy, 
where = (— a,0, 0, 0) [and = g^^^n^ = (1/a, —(5'^/a)\ is the unit normal vector of the spacelike 
hypersurface, E (see Figure 2). By introducing the extrinsic curvature, 

Kij = ~£nl^J, (2.2) 

and using the Gauss-Codacci relation, the Hamiltonian density of the Einstein equations can be written 
as 

Hgr = 7r% -A where C = ^/^R = a^[^^^R - + KijK% (2.3) 
where tt*^ is the canonically conjugate momentum to jij, 

^ij = IL = -^{K'i - Kf^), (2.4) 

O'Jij 

omitting the boundary terms. The variation of ficR with respect to a and Pi yields the constraints, 
and the dynamical equations are given by jij = and tt*-^ = — ^j^^- 



Figure 2: Concept of time evolution of space-time: foliations of 3-dimensional hypersurface. The lapse 
and shift functions are often denoted q or A^, and /3* or AT*, respectively. (This figure is missing in 
gr-qc version due to the limitation of the file size.) 
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2.0.2 The standard ADM formulation 

In the version of Smarr and York, Kij was used as a fundamental variable instead of the conjugate 
momentum tt*-' (see also the footnote ^). 



The Standard ADM formulation [80, 98]: Box 2.1 

The fundamental dynamical variables are {'^ij^Kij), the three-metric and extrinsic curvature. 
The three-hypersurface E is foliated with gauge functions, (a,/3*), the lapse and shift vector. 

• The evolution equations: 

dt-iij = -2aKij + Dipj + Djf3„ (2.5) 
dtKij = a ^^'^Rij + aKKij - 2aKikK^ j - DiDja 

+{Di(5'')Kkj + {D,f]'')Kki + fi^D^Kij (2.6) 

where K = K\, and ^^^Rij and Di denote three-dimensional Ricci curvature, and a 
covariant derivative on the three-surface, respectively. 

• Constraint equations: 

y_ADM _ i^)R + K^ -KijK'^ ^Q, (2.7) 
^ADM ._ DjK^i- DiK kO, (2.8) 

where ^^^R =^^^ i?\: these are called the Hamiltonian (or energy) and momentum con- 
straint equations, respectively. 



The formulation has 12 free first-order dynamical variables (jijjKij), with 4 freedom of gauge choice 
{a, Pi) and with 4 constraint equations, (2.7) and (2.8). The rest freedom expresses 2 modes of 
gravitational waves. 

The constraint propagation equations, which are the time evolution equations of the Hamiltonian 
constraint (2.7) and the momentum constraints (2.8), can be written as 



The Constraint Propagations of the Standard ADM: 


Box 2.2 


dtH = P^{djH) + 2aKH-2af^{diMj) 




+a(5amfe)(27"^S'^' - ^""'I'nMj - 4Y^{dja)Mi, 


(2.9) 


dtMi = -{l/2)a{d{H)-{dia)n + P\djMi) 




+aKMi - /3''7'\daik)Mj + {dAh''' Mj. 


(2.10) 


Further expressions of these constraint propagations are in Appendix A. 





Prom these equations, we know that if the constraints are satisfied on the initial slice S, then the 
constraints are satisfied throughout evolution. The normal numerical scheme is to solve the elliptic 

^We remark that there is one replacement in (2.6) using (2.7) in the process of conversion from the original ADM to 
the standard ADM equations. This is the key issue in the later discussion, and we shall be back this point in §3.2. 
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constraints for preparing the initial data, and to apply the free evolution (solving only the evolution 
equations). The constraints are used to monitor the accuracy of simulations. 



2.0.3 Remarks 

The ADM formulation was the standard formulation for numerical relativity up to the middle 90s. 
Numerous successful simulations were obtained for the problems of gravitational collapse, critical 
behavior, cosmology, and so on. However, stability problems have arisen for the simulations such 
as the gravitational radiation from compact binary coalescence, because the models require quite a 
long-term time evolution. 

The origin of the problem was that the above statement in Italics is true in principle, but is not 
always true in numerical applications. A long history of trial and error began in the early 90s. From 
the next subsection we shall look back of them by summarizing "three strategies". We then unify 
these three roads as "adjusted systems", and as its by-product we show in §3.2 that the standard 
ADM equations has a constraint violating mode in its constraint propagation equations even for a 
single black-hole (Schwarzschild) spacetime [78]. 

2.1 Strategy 1: Modified ADM formulation by Nakamura et al (The BSSN for- 
mulation) 

2.1.1 Introduction 

Up to now, the most widely used formulation for large scale numerical simulations is a modified ADM 
system, which is now often cited as the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation. 
This reformulation was first introduced by Nakamura et al. [63, 64, 74]. The usefulness of this reformu- 
lation was re-introduced by Baumgarte and Shapiro [15], then was confirmed by another group to show 
a long-term stable numerical evolution [5, 7]. The procedure is to apply conformally decomposition 
of the ADM variables and to implement their dynamical equations with several replacements. 

2.1.2 Basic variables and equations 

The widely used notation[15] introduces the variables {ip,^ij,K,Aij,r'^) instead of {-jijjKij), where 

<p = (l/12)log(dct7,,), (2.11) 

lij = e-^^jij, (2.12) 

K = f^Kij, (2.13) 

A,j = e-^'^{K,j-{l/3)jijK), (2.14) 

r = tik^"- (2.15) 

The new variable was introduced in order to calculate Ricci curvature more accurately. F* also 
contributes to making system re-produce wave equations in its linear limit. 
In BSSN formulation, Ricci curvature is not calculated as 



but 



Rfj^'^^ - dkTfj - d,Tij + F-jFffc - F^^^Ffj, (2.16) 



Rl = -2DiDjip-2%D''Dkcp + A{Di^){Dj^)-A%{D''cp){Dkcp), (2.18) 

■ l{i 



Rij = -(l/2)7''=ai5fe7i,-+%(i5,)f'= + f'=f(,,)fc + 27'-ff(.f,)^^ (2.19) 
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where Di is covariant derivative associated with ^jj. These are approximately equivalent, but Rf^^^^ 
does have wave operator apparently in the flat background limit, so that we can expect more natural 
wave propagation behavior. 

Additionally, the BSSN requires us to impose the conformal factor as 

7(:= det7,j) = 1, (2.20) 

during evolution. This is a kind of definition, but can also be treated as a constraint. 



The BSSN formulation [63, 64, 74, 15]: Box 2.3 

The fundamental dynamical variables are ((^, 7ij,K,JLjj,P). 

The three-hypersurface S is foliated with gauge functions, (a,/3*), the lapse and shift vector. 

• The evolution equations: 

= -{l/Q)aK + {l/Q)P\di^) + {d^p'), (2.21) 

df^ij = -2«4- + 7ifc(5i/3') + ljk{dil5^) - (2/3)7i,(5fe/3'=) + /?'(5fc7ii), (2.22) 

dfK = -DWia + aAijA'^ + {l/3)aK^ +p\diK), (2.23) 
5^4- = -e-^^(A^,-a)^^ + e-^^a(i?f^^)^^ + aif4--2ai,fci^- 

+{^^P'')AkJ + i^,p')Ak^ - (2/3) (9^/3*^)4- + P'idkA^), (2.24) 
dfr = -2{dja)A'^ + 2a(f;fe#^' - {2/3)f^{djK) + 6A'^{djip)) 

-djiP'^idkf^) - i^Kdkd') - f\dkl3^) + {2/3)r{dkP')). (2.25) 

• Constraint equations: 



q^BSSN 




(2.26) 


j^BSSN 


= Mf"^, 


(2.27) 




= t'-f'fi„ 


(2.28) 


A 


= Ajf', 


(2.29) 


s 


= 7-1. 


(2.30) 



(2.26) and (2.27) are the Hamiltonian and momentum constraints (the "kinematic" con- 
straints), while the latter three are "algebraic" constraints due to the requirements of 
BSSN formulation. 



Hereafter we will write T-l^^^^ and M^^^^ simply as H. and M respectively. 

Taking careful account of these constraints, (2.26) and (2.27) can be expressed directly as 

rC = e-^'PR-8e-^'PDWj<f-8e-*'^{D^(p){Dj(p) + {2/3)K^-AijA'^-{2/3)AK, (2.31) 
Mi = QAh{D,^)-2A{Di^)-{2/3){D^K) + ^^^{b,Ak^. (2.32) 

In summary, the fundamental dynamical variables in BSSN are {(p,^ij, i^,Ajj,r*), 17 in all. The 
gauge quantities are (q;,/3*) of which there are 4, and the constraints are {H,Aii,G^,A,S), i.e. 9 
components. As a result, 4 (2 by 2) components are left which correspond to two gravitational 
polarization modes. 
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2.1.3 Remarks 



Why is BSSN better than the standard ADM? Together with numerical comparisons with the standard 
ADM case[7, 57], this question has been studied by many groups using different approaches. Using 
numerical test evolution, Alcubierre et al [5] found that the essential improvement is in the process of 
replacing terms by the momentum constraints. They also pointed out that the eigenvalues of BSSN 
evolution equations have fewer "zero eigenvalues" than those of ADM, and they conjectured that the 
instability might be caused by these "zero eigenvalues" . Miller [60] applied von Neumann's stability 
analysis to the plane wave propagation, and reported that BSSN has a wider range of parameters, 
which produces stable evolution. Analogical conformal decomposition of the Maxwell equations are 
also reported [53]. An effort was made to understand the advantage of BSSN from the point of 
hyperbolization of the equations in its linearized limit [5, 68]. These studies provide some support 
regarding the advantage of BSSN, while it is also shown an example of an ill-posed solution in BSSN 
(as well in ADM) by Prittelli and Gomez [40]. (Inspired by the BSSN's conformal decomposition, 
several related hyperbolic formulations have also been proposed [11, 42, 8].) 

As we shall discuss in §3.3, the stability of the BSSN formulation is due not only to the introductions 
of new variables, but also to the replacement of terms in the evolution equations using the constraints. 
Further, we will show several additional adjustments to the BSSN equations which are expected to 
give us more stable numerical simulations. 

Recently, Laguna and Shoemaker [55] modified the BSSN slightly, and found a great improvement 
in simulating a Schwarzschild black-hole. 



The Laguna-Shoemaker version of BSSN [55]: Box 2.4 

Modifications to BSSN are 

• to introduce conformal scalings also to K, Aij and N as 

where jij = e~^'^ jij and n is a parameter (n = recovers the BSSN variables) 

• to use a mixed indices form ^^^^A^j rather than A^j. 

• to use a densitized lapse N rather than a. 



There is no explicit explanation why these small changes work better than before, but we expect that 
our method can also be applied to finding the reason. 

2.2 Strategy 2: Hyperbolic reformulations 

2.2.1 Definitions, properties, mathematical backgrounds 

The second effort to re-formulate the Einstein equations is to make the evolution equations reveal 
a first-order hyperbolic form explicitly. This is motivated by the expectation that the symmetric 
hyperbolic system has wcll-posed properties in its Cauchy treatment in many systems and also that 
the boundary treatment can be improved if we know the characteristic speed of the system. As a 
comprehensive review of the hyperbolic formulation, we refer to those by Choquet-Bruhat and York 
(1980) [31], Geroch (1996) [43], Reula (1998) [67], and Priedrich-Rendah (2000) [37], 
We use the following definition: 
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Hyperbolic formulations Box 2.5 

We say that the system is a first-order (quasi-linear) partial differential equation system, if a 
certain set of (complex- valued) variables Mq (a = 1, • • • , n) forms 

dtUa = M^^a{u) diU0+Ma{u), (2.33) 

where M (the characteristic matrix) and M are functions of u but do not include any derivatives 
of u. Further we say the system is 

• a weakly hyperbolic system, if all the eigenvalues of the characteristic matrix are real. 

• a strongly hyperbolic system (or a diagonalizable / symmetrizable hyperbolic system), if 
the characteristic matrix is diagonalizable (has a complete set of eigenvectors) and has all 
real eigenvalues. 

• a symmetric hyperbolic system, if the characteristic matrix is a Hermitian matrix. 



We treat A^'^^ as a n x n matrix when the i-index is fixed. The following properties of these matrices 
apply to every basis of Z-index. We say A' is an eigenvalue of M^^a when the characteristic equation, 
det(7W'^a - A'J'^a) = 0, is satisfied. The eigenvectors, p°', are given by solving = X^p^°'. 

The strong hyperbolicity is identified, e.g. by judging whether the multiplicity of its eigenvalue, nx, 
satisfies rank(J^^^ ^ — ^^^^a) = n — n\ for all A. In order to define the symmetric hyperbolic system, 
we need to declare an inner product {u\u) to determine whether M''^a is Hermitian. In other words, 
we are required to define the way of lowering the index a of u". We say M^^a is Hermitian with 
respect to this index rule, when A^'/3q, = M^ap for every I, where the overhead bar denotes complex 
conjugate. To avoid this requirement of the definition of the inner product, people sometimes use 
the word symmetrizable, when the characteristic matrix becomes Hermitian by a certain symmetrizer 
(positive definite matrix). In our classification, this is only equivalent to the strongly hyperbolic 
system. ^ 

Writing the system in a hyperbolic form is a quite useful step in proving that the system is well- 
posed. The mathematical well-posedness of the system means (1°) local existence (of at least one 
solution u), (2°) uniqueness (i.e., at most solutions), and (3°) stability (or continuous dependence of 
solutions {u} on the Cauchy data) of the solutions. The resultant statement expresses the existence 
of the energy inequality on its norm, 

I \u{t) 1 1 < e'"'' I |u(t = 0) 1 1 , where < t < t, a = const. (2.34) 

This indicates that the norm of u{t) is bounded by a certain function and the initial norm. Remark 
that this mathematical does not mean that the norm u{t) decreases along the time evolution. 
The inclusion relation of the hyperbolicities is, 

symmetric hyperbolic C strongly hyperbolic C weakly hyperbolic. (2.35) 

The Cauchy problem under weak hyperbolicity is not, in general, C°° well-posed. At the strongly 

hyperbolic level, we can prove the finiteness of the energy norm if the characteristic matrix is indepen- 
dent of u (cf [81]), that is one step definitely advanced over a weakly hyperbolic form. Similarly, the 

Several groups use a slightly different definition for a symmetric hyperbolic system: defining when the principal 
symbol of the system iM''^ aki is anti-Hermitian for arbitrary vector ki. The two definitions are equivalent when the 
vector fc; is real-valued, but different eigenvalues. In our definition, all eigenvalues are real-valued, while in the other 
they are all pure imaginary. 
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well-posedness of the symmetric hyperbolic is guaranteed if the characteristic matrix is independent of 
u, while if it depends on u we have only limited proofs for the well-posedness. From the mathematical 
point of view, proving well-posedness with less strict conditions is an old but still active research 
problem. Therefore we have to recall that even if we construct a symmetric hyperbolic system in 
general relativity, that equation does not necessarily guarantee numerical stability. 

From the point of numerical applications, to hyperbolize the evolution equations is quite attractive, 
not only for its mathematically well-posed features. The expected additional advantages are the 
following. 

(a) It is well known that a certain flux conservative hyperbolic system is taken as an essential formu- 
lation in the computational Newtonian hydrodynamics when we control shock wave formations 
due to matter (e.g. [45]). 

(b) The characteristic speed (eigenvalues of the principal matrix) is supposed to be the propagation 
speed of the information in that system. Therefore it is naturally imagined that these magnitudes 
are equivalent to the physical information speed of the model to be simulated. 

(c) The existence of the characteristic speed of the system is expected to give us an improved 
treatment of the numerical boundary, and/or to give us a new well-defined Cauchy problem 
within a finite region (the so-called initial boundary value problem, IBVP). 

These statements sound reasonable, but have not yet been generally confirmed in actual numerical 
simulations. But we are safe in saying that the formulations are not yet well developed to test these 
issues. For example, IBVP studies are preliminary yet, and most works are based on a particular 
symmetric hyperbolic and in a limited space-time symmetry [26, 29, 38, 49, 81, 82, 83, 84]. We will 
come back to this issue in §2.2.3 or §4, but meanwhile let us view the hyperbolic formulations from 
the comparisons of pure evolution equations. 

We note that rich mathematical theories on partial differencial equations are obtained in a first- 
order form, while there is a study on linearized ADM equations in a second-order form [54]. 

2.2.2 Hyperbolic formulations of the Einstein equations 

As was discussed by Geroch [43], most physical systems can be expressed as symmetric hyperbolic 
systems. In order to prove that the Einstein's theory is a well-posed system, to hyperbolize the 
Einstein equations is a long-standing research area in mathematical relativity. As we mentioned in 
the introduction, numerical relativity shed light on this mathematical problem. 

The standard ADM system does not form a first order hyperbolic system. This can be seen 
immediately from the fact that the ADM evolution equation (2.6) has Ricci curvature in RHS. So far, 
several first order hyperbolic systems of the Einstein equation have been proposed. In constructing 
hyperbolic systems, the essential procedures are (1°) to introduce new variables, normally the spatially 
derivatived metric, (2°) to adjust equations using constraints. Occasionally, (3°) to restrict the gauge 
conditions, and/or (4°) to rescale some variables. Due to process (1°), the number of fundamental 
dynamical variables is always larger than that of ADM. 

In the following discussion, we briefly review several hyperbolic systems of the Einstein equa- 
tions. We only mention the systems which applied numerical comparisons. See Table. 1 for a more 
comprehensive list. 

• The Bona-Masso formulation [18, 19, 20, 21] 

This was the first active effort to apply hyperbolized equations to numerical relativity. They 

introduced auxiliary variables to reduce the system first order in space: = dklna, = 
(l/2)9fc/3*, and Dfcjj = {l/2)dkgij, and construct a first order flux conservative system. In their 
latest formulation [21], the set of dynamical variables is (aj'jij, Kij, Ai, Drij,Vi = 2D^irf), 37 
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functions in total, and the lapse function is restricted to a certain functional condition. The 
system is a symmetrizable hyperbolic. They observed improved numerical stability compared to 
that of the standard ADM system in spherical symmetric spacetime evolution [20] . An advantage 
of having characteristic speed is also applied to improve the treatment of the outer boundary 
condition [22]. However, the appearance of shock formation was also reported unless the lapse 
function is determined by solving the elliptic or parabolic equations [2, 4]. 

The Einstein- Ricci system [32, 1] / Einstein-Bianchi system [9] 

Series of works by Choquet-Bruhat, York, and their colleagues developed hyperbolic systems in 
slightly different ways. They intended to construct wave-type equations (D^ = S) for (jijjKij), 
that require the introduction of new variables (45 or 66 variables in total) and the use of Bianchi 
identities. The resultant system is a third-order system in time for since they use the equation 
dtRij, but has only physical characteristic speed (zero or light speed). Scheel et al [70, 71] devel- 
oped a numerial code with this formulation, and reported that they could evolve Schwarzschild 
black hole (with 1-dimcnsional code) quite long time (~ lO^M), but futher modifications to 
equations were ncccsarrily. 

The Einstcin-Christoffel system [10] 

Anderson and York [10] constructed a symmetrizable hyperbolic system which only has physical 
characteristic speed. The variables are fundamentally {jij, Kij,Tkij), and they also derived a set 
of equations with {'jij, Kij,Gkij = dkjij), 30 functions. Their formulation differs from the Bona- 
Masso formulation in the way the momentum constraint is used to make the system hyperbolic. 
Using a model of plane wave propagation, Bardeen and Buchman [17] numerically compared this 
formulation with ADM and Bona-Masso with slight changes in variables. Their experiments are 
of plane-symmetric wave propagations, and they studies the boundary treatment in detail. We 
will mention their work later. 

The Ashtekar formulation [13] 

Ashtekar's reformulation of space-time was introduced to provide a new non-perturbative ap- 
proach to quantum gravity. The new basic variables are the densitized inverse triad, and 
the S0(3,C) self-dual connection, Af, where the indices indicate the 3-spacetime, and 

a,b,--- are for SO (3) space. The formulation requires additional constraints and reality con- 
ditions in order to describe the classical Lorentzian space-time evolution, but the evolution 
system itself forms a weakly hyperbolic system. By keeping the number of variables the same 
[(£^*,^")=total 18 (minimum ever)], we can construct both strongly and symmetric hyperbolic 
systems by adjusting equations with constraints and/or restricting gauge conditions [92, 93]. 
The authors made numerical comparisons between the hyperbolicity of the systems with plane 
symmetric gravitational wave propagation using a periodic boundary condition [77, 94]. The 
outcome is that there are no drastic differences in numerical stability between the three levels 
of the hyperbolic equations. We will describe the details in Appendix B. 

The Prittelli-Reula formulation [41, 81] 

This is a one-parameter family of a symmetric hyperbolic system. The procedure is to define 
new variables, to adjust evolution equations with constraints (with one parameter), and to 
densitize the lapse function (with one parameter) . The variables are [jij , M*-'^ = ^ (7*-' + 
dl^-^lrsY^ k)^ = + bY^K) where a,b are two other parameters, and make a total of 30. 
They define a symmetric hyperbolic system by non-diverging energy norm, and restricted free 
parameter spaces. A numerical comparison was also made for Gowdy spacetime [44], with quite 
similar conclusions as in the Ashtekar version. 

The Conformal Field equations [36] 

A series of works by Friedrich [36] attempted to construct a 3-|- 1 formulation with hyperboloidal 
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foliations (i.e. asymptotically null foliations), and with comformal compactification. This is the 
ultimate plan to remove the outer boundary problem in numerical simulation, and to provide a 
suitable foliation for gravitational radiation problem. However, the current equations are rather 
quite complicated. In its metric-based expression [46], the evolution variables are 57; -yijjKij, 
the connection coefficients j"'bc, projections ^^'^^Ra = n^^a'^Rhc and ^^'^^Rah = Ja^Jb'^Rbd of 4- 
dimensional Ricci tensor Rab, the electric and magnetic components of the rescaled Weyl tensor 
Cabc^i and the comformal factor $7 and its related quantities = n"Vai7, Va^^, V"Vai7. By 
specifying suitable gauge functions (a, R) where R is the Ricci scalar, then the total system 
forms a symmetric hyperbolic system. Applications to numerical relativity are in progress, but 
have not yet reached the stage of applying evolution in a non-trivial metric. For more details, 
see reviews e.g. by Frauendiener [34] or by Husa [47]. 

• The Kidder-Scheel-Teukolsky (KST) formulation [52] 

This set of equations is a generalized formulation of the previous ones. It has 12 free parameters, 
and includes both Anderson- York [10] and Frittelli-Reula[41] formulations as a subset. Therefore 
it might be useful to discuss further hyperbolization methods starting with this KST formulation. 
We briefly summarize it in Box 2.6. 



The Kidder-Scheel-Teukolsky (KST) formulation [52]: Box 2.6 

• Starting from a set of {gij , Kij , dkij = OkQij), the fundamental dynamical variables are 
defined {gij , Pij , M^ij) , as 

P,j = K,j + zgijK, (2.36) 
Mkij = {l/2)[kdkij + ed(ij)k + gij{ddk + hhk) + gk(i{cdj) + dhj))], (2.37) 

where dk = g°'^dkab-> bk = g°'^dabk-> a^d {d,b,c,d,e, k,z) are "kinematical" parameters. 

• The 3-hypersurface E is fohated with gauge functions, (a,/?*), the lapse and shift vector. 
However, the densitized lapse, Q = log{ag~'^), (with a parameter a) is actually used. 

• The evolution equations are adjusted with constraints [in a version of {gij , Kij , dkij)] 

-2aKij, (2.38) 
(•••) + lag.jH + Ca5"'C„(,,)fe, (2.39) 
(•••) + mgk{iMj) + xagijMk, (2.40) 

where do = dt — £[3 and (7,C)^)X) are parameters. The terms (■ ■ ■) arc original terms 
derived from the ADM equations and are available as (2.14) and (2.24) in [52]. 

• Constraints are {H,Mi,Ckiij), where CkUj = d\kdi-\ij. 



In short, the number of dynamical variables and constraints are 30 and 22 and there are 12 free 
parameters. Although there is no specific method to specify 12 free parameters KST showed numerical 
examples of the Schwarzschild black hole (mass M) evolution, which runs quite a longer time {t ~ 
6000M) [52, 72]. 



dogij = 

doKij = 
dodkij = 
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We think the essential advantage of the KST system is the introduction of "kinematical" parame- 
ters. These 6 parameters 

• do not change the eigenvalues of evolution eqs., 

• do not affect the principal part of the constraint evolution eqs., 

• do affect the eigenvectors of the evolution system, and 

• do affect the nonlinear terms of evolution eqs/constraint evolution eqs. 

As Calabrese et al [27] pointed out, KST equations at linearized level on the fiat spacetime have no 
non-principal terms, and these "kinematical" parameters finally enable ns to discuss the features of 
hyperbolicity in numerical experiments. Several variations of the KST formulation are presented in 
[69]. Recently, Lindblom and Scheel [58] tried to explain the relation between the numerical error 
growing rate and the predicted error growing rate from the characteristic matrix of KST evolution 
equations. Their trial did not succeed in matching the two completely, but at least began to reveal a 
similar order before non-linear blow-up begins. 

2.2.3 Remarks 

When we discuss hyperbolic systems in the context of numerical stability, the following questions 
should be considered: 

Questions to hyperbolic formulations on its applications to numerics 

Q(A) Prom the point of the set of evolution equations, does hyperbolization actually contribute 
to numerical accuracy and stability? Under what conditions/situations will the advantages 
of hyperbolic formulation be observed? 

Q(B) If the answer to Q(A) is affirmative, which level of hyperbolicity is practically useful for 
numerical applications? Strongly hyperbolic, symmetric hyperbolic, or other? 

Q(C) If the answer to Q(A) is negative, then can we find other practical advantages to hyper- 
bolization? Treatment of boundary conditions, or other? 

Unfortunately, we do not have conclusive answers to these questions, but many experiences are 
being accumulated. Several earlier numerical comparisons reported the stability of hyperbolic formu- 
lations [20, 21, 22, 70, 71]. But we have to remember that this statement went against the standard 
ADM formulation, which has a constraint-violating mode for Schwarzschild spacetime as has been 
shown recently (see §3.2). 

These partial numerical successes encouraged the community to formulate various hyperbolic sys- 
tems. Recently, Calabrese et al [28] reported there is a certain differences in the long-term convergence 
features between weakly and strongly hyperbolic systems on the Minkowskii background space-time. 
However, several numerical experiments also indicate that this direction is not a complete success. 

Objections from numerical experiments 

• Above earlier numerical successes were also terminated with blow-ups. 

• If the gauge functions are evolved according to the hyperbohc equations, then their finite 
propagation speeds may cause pathological shock formations in simulations [2, 4]. 

• There are no drastic differences in the evolution properties between hyperbolic systems 

(weakly, strongly and symmetric hyperbolicity) by systematic numerical studies by Hern 
[44] based on Frittelli-Reula formulation [41], and by the authors [77] based on Ashtekar's 
formulation [13, 92, 93]. 

• Proposed symmetric hyperbolic systems were not always the best ones for numerical evo- 
lution. People are normally still required to reformulate them for suitable evolution. Such 
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efforts are seen in the applications of the Einstein- Ricci system [71], the Einstein-Christoffel 
system [17], the conformal field equations [47], and so on. 

• Bardeen and Buchmann [17] confirmed the usefulness of hyperbolicity when they treated 
numerical boundary conditions in plane-symmetric wave propagation problem, but they 
also mentioned that the same techniques can not be applied in 2 or 3-dimensional cases. 

Of course, these statements only casted on a particular formulation, and therefore we have to be 
careful not to over-emphasize the results. In order to figure out the reasons for the above objections, 
it is worth stating the following cautions: 
Remarks on hyperbolic formulations 

(a) Rigorous mathematical proofs of well-posedness of PDE are mostly for simple symmetric or 
strongly hyperbolic systems. If the matrix components or coefficients depend on dynamical 
variables (as in all any versions of hyperbolized Einstein equations), almost nothing was 
proved in more general situations. 

(b) The statement of "stability" in the discussion of well-posedness refers to the bounded 
growth of the norm, and does not indicate a decay of the norm in time evolution. 

(c) The discussion of hyperbolicity only uses the characteristic part of the evolution equations, 
and ignores the rest. 

We think the origin of confusion in the community results from over-expectation on the above 
issues. Mostly, point (c) is the biggest problem, as was already pointed out in several places [94, 56, 78]. 
The above numerical claims from Ashtekar and Frittelli-Reula formulations were mostly due to the 
contribution (or interposition) of non-principal parts in evolution. 

Regarding this issue, the recent KST formulation finally opens the door. As we saw, KST's 
"kinematic" parameters enable us to reduce the non-principal part, so that numerical experiments 
are hopefully expected to represent predicted evolution features from PDE theories. At this moment, 
the agreement between numerical behavior and theoretical prediction is not yet perfect but close [58]. 
If further studies reveal the direct correspondences between theories and numerical results, then the 
direction of hyperbolization will remain as the essential approach in numerical relativity, and the 
related IBVP researches will become a main research subject in the future. 

Meanwhile, it will be useful if we have an alternative procedure to predict stability including the 
effects of the non-principal parts of the equations, which are neglected in the discussion of hyperbolicity. 
Our proposal of adjusted system in the next subsection may be one of them. 

2.3 Strategy 3: Asymptotically constrained systems 

The third strategy is to construct a robust system against the violation of constraints, such that the 
constraint surface is an attractor (Figure 3). The idea was first proposed as "A-system" by Brodbeck 
et al [24], and then developed in more general situations as "adjusted system" by the authors [94]. 

2.3.1 The "A-system" 

Brodbeck et al [24] proposed a system which has additional variables A that obey artificial dissipative 
equations. The variable As are supposed to indicate the violation of constraints and the target of the 
system is to get A = as its attractor. Their proposal can be summarized as in Box 2.7. 
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The "A-system" (Brodbeck-Frittelli-Hiibner-Reula) [24]: Box 2.7 

For a symmetric hyperbolic system, add additional variables A and artificial force to reduce the 
violation of constraints. 
The procedure: 



1. 
2. 

3. 
4. 



Prepare a symmetric hyperbolic evolution system dtu = MdiU + 



Introduce A as an indicator of violation of 
constraint which obeys dissipative eqs. of motion 

Take a set of (n, A) as dynamical variables 

Modify evolution eqs so as to form 
a symmetric hyperbolic system 



dt 
dt 



dtX = aC-pX 
(q / 0,/3 > 0) 

A/ ^ I F 



A F 
F 



The main ideas are to introduce additional variables, As, to impose dissipative dynamical equations for 
As, and to construct a symmetric hyperbolic system for both original variables and As. Since the total 
system is designed to have symmetric hyperbolicity, the evolution is supposed to be unique. Brodbeck 
et al showed analytically that such a decay of As can be seen for sufficiently small A(> 0) with a choice 
of appropriate combinations of as and /3s. 

Brodbeck et al presented a set of equations based on Frittelli-Reula's symmetric hyperbolic for- 
mulation [41]. The version of Ashtekar's variables was presented by the authors [76] for controlling 
the constraints or reality conditions or both (see §B.2.2). The numerical tests of both the Maxwell- 
A-system and the Ashtekar- A-system were performed [94], and confirmed to work as expected (see 
§B.3.3). Although it is questionable whether the recovered solution is true evolution or not [79], we 
think the idea is quite attractive. To enforce the decay of errors in its initial perturbative stage seems 
the key to the next improvements, which are also developed in the next section on "adjusted systems" . 

However, there is a high price to pay for constructing a A-system. The A-system can not be 
introduced generally, because (i) the construction of A-system requires the original evolution equations 
to have a symmetric hyperbolic form, which is quite restrictive for the Einstein equations, (ii) the final 
system requires many additional variables and we also need to evaluate all the constraint equations at 
every time step, which is a hard task in computation. Moreover, (iii) it is not clear that the A-system 
is robust enough for non-linear violation of constraints, or that A-system can control constraints which 
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Constrained/ Surface 
(satisfies /Einstein's constraints) 




time 



Figure 3: Schematic picture of "asymptotically constrained system". 
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do not have any spatial differential terms. 
2.3.2 The "adjusted system" 

Next, we propose an alternative system which also tries to control the violation of constraint equations 
actively, which we named "adjusted system". We think that this system is more practical and robust 
than the previous A-system. 



The Adjusted 


system (essentials) [94]: Box 2.8 


Purpose: 


Control the violation of constraints by reformulating the system so as to 
have a constrained surface as attr actor. 


Procedure: 


Add a particular combination of constraints to the evolution equations, 
and adjust its multipliers. 


Theoretical support: Eigenvalue analysis of the constraint propagation equations. 


Advantages: 


Available even if the base system is not symmetric hyperbolic. 


Advantages: 


Keeps the number of the variables the same as in the original system. 



We will describe the details in the next section, but summarize the procedures in advance: 



The Adjusted system (procedures): 




Box 2.9 


1. Prepare a set of evolution eqs. 


dtu = 


JdiU + K 


2. Add constraints in RHS 


dtu = 


JdiU + K +kC 


Choose the coeff. k so as to make the 
3. eigenvalues of the homogenized adjusted 
dfC eqs negative reals or pure imaginary. 
(Sec Box 3.2 and 3.3) 


dtC 
dtC 


= Dd,,C + EC 

= DdiC + EC +FdiC + GC 
^ V ' 


The details are in §3. 







The process of adjusting equations is a common technique in other re-formulating efforts as we re- 
viewed. However, we try to employ the evaluation process of constraint amplification factors as an 
alternative guideline to hyperbolization of the system. 

For the Maxwell equation and the Ashtekar version of the Einstein equations, we numerically 
found that this idea works to reduce the violation of constraints, and that the effects are much better 
than by constructing its symmetric hyperbolic versions [94] (see §B.3.4). The idea was apphed to the 
standard ADM formulation which is not hyperbolic and several attractive adjustments were proposed 
[95, 78] (see §3.2). This analysis was also applied to explain the advantages of the BSSN formulation, 
and again several alternative adjustments to BSSN equations were proposed [96] (see §3.3). We will 
explain these issues in the next section. 
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3 A unified treatment: Adjusted System 



This section is devoted to present our idea of "asymptotically constrained system" , which was briefly 
introduced in Box 2.8 and 2.9 in the previous section. We begin with an overview of the adjusting 
procedure and the idea of background structure in §3.1. Next, we show the applications both to ADM 
(§3.2) and BSSN (§3.3) formulations. Original references can be found in [94, 95, 78, 96]. 

3.1 Procedures : Constraint propagation equations and Proposals 

Suppose we have a set of dynamical variables u°'{x^,t), and their evolution equations, 

att." = /(n«,5in«,---), (3.1) 

and the (first class) constraints, 

C°^(u",aiti",---) ~ 0. (3.2) 

Note that we do not require (3.1) forms a first order hyperbolic form. We propose to investigate the 
evolution equation of C"^ (constraint propagation), 

5tC° = 5(C°,5iC",---), (3.3) 

for predicting the violation behavior of constraints in time evolution. We do not mean to integrate 
(3.3) numerically together with the original evolution equations (3.1), but mean to evaluate them 
analytically in advance in order to reformulate the equations (3.1). 

There may be two major analyses of (3.3); (a) the hyperbolicity of (3.3) when (3.3) is a first order 
system, and (b) the eigenvalue analysis of the whole RHS in (3.3) after a suitable homogenization. 

(a) Hyperbolicity analysis of (3.3) 

If (3.3) forms a first-order system, the standard hyperbolic PDE analysis is applicable. As we 
viewed in §2.2, the analysis is mainly to identify the level of hyperbolicity and to calculate the 
characteristic speed of the system, from eigenvalues of the principal matrix. 

For example, the evolution equations of the ADM formulations, (3.1) [(2.5) and (2.6)], do not 
form a first-order system, while their constraint propagation equations, (3.3) [(2.9) and (2.10)], 
do form. Therefore, we can apply the classification on the hyperbolicity (weakly, strongly or 
symmetric) to the constraint propagation equations. However, if one adjusts the ADM equations 
with constraints, then this first-order characters will not be guaranteed. 

As we mentioned in §2.2.3, another big problem in the hyperbolic analysis is that it only discusses 
the principal part of the system. If there is a method to characterize the non-principal part, 
then that will help to clarify our understanding of evolution behavior. 



(b) Eigenvalue analysis of the whole RHS in (3.3) after a suitable homogenizing procedure. 



This analysis may compensate for the above problems. 


Amplification Factors of Constraint Propagation equations: 

We propose to homogenize (3.3) by a Fourier transformation, e.g. 


Box 3.1 






where C(x, tf = j C{k, tf exp{ik ■ x)d^k, 


(3.4) 


then to analyze the set of eigenvalues, say As, of the coefficient matrix, M"^, 
call As the constraint amplification factors (CAFs) of (3.3). 


in (3.4). We 
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The CAFs predict the evolution of constraint violations. We therefore can discuss the "dis- 
tance" to the constraint surface using the "norm" or "compactness" of the constraint violations 
(although we do not have exact definitions of these "• • •" words). 

The next conjecture seems to be quite useful to predict the evolution feature of constraints: 



Conjecture on Constraint Amplification Factors (CAFs): Box 3.2 

(A) If CAF has a negative real-part (the constraints are forced to be diminished), then 
we see more stable evolution than a system which has positive CAF. 

(B) If CAF has a non-zero imaginary-part (the constraints are propagating away), then 
we see more stable evolution than a system which has zero CAF. 



We found that the system becomes more stable when more As satisfy the above criteria. (The 
first observations were in the Maxwell and Ashtekar formulations [77, 94]). Actually, supporting 
mathematical proofs are available when we classify the fate of constraint propagations as follows. 



Classification of Constraint propagations: Box 3.3 

If we assume that avoiding the divergence of constraint norm is related to the numerical 
stability, the next classification would be useful: 

(CI) Asymptotically constrained : All the constraints decay and converge to zero. 
This case can be obtained if and only if all the real parts of CAFs are negative. 

(C2) Asymptotically hounded : All the constraints are bounded at a certain value, (this 
includes the above asymptotically constrained case.) 

This case can be obtained if and only if (a) all the real parts of CAFs are not positive 
and the constraint propagation matrix M"^ is diagonalizable, or (b) all the real parts 
of CAFs are not positive and the real part of the degenerated CAFs is not zero. 

(C3) Diverge : At least one constraint will diverge. 

The details are shown in [97] . 



This classification roughly indicates the heuristic statements (A) in Box 3.2. A practical proce- 
dure for this classification is drawn in Figure 4. 

We remark that this eigenvalue analysis requires the fixing of a particular background spacetime, 
since the CAFs depend on the dynamical variables, u". 

The above features of the constraint propagation, (3.3), will differ when we modify the original 
evolution equations. Suppose we add (adjust) the evolution equations using constraints 

dtu'' = /(u", diu'', • • •) + FiC, diC'', ■ ■ •), (3.5) 

then (3.3) will also be modified as 

^tC" = 5(C", a^C", • • •) + G(C", • • •). (3.6) 

Therefore, the problem is how to adjust the evolution equations so that their constraint propagations 
satisfy the above criteria as much as possible. 
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Ql: Is there a CAF which real part is positive? 

NO / YES ► Diverge 

I 

Q2: Are all the real parts of CAFs negative? 

Asymptotically 
' " Contained 

I 

Q3: Is the constraint propagation matrix diagonalizable? 

NO / YES ► Asymptotically 

^ Bounded 

Q4: Is a real part of the degenerated CAFs is zero? 

"^^S / NO ► Asymptotically 

^ Bounded 

Q5: Is the associated Jordan matrix diagonal? 

NO / YES ► Asymptotically 

^ Boimded 

Diverge 

Figure 4: A flowchart to classify the constraint propagations (Box 3.3). 
3.2 Adjusted ADM formulations 

We show an application of the above idea, to the standard ADM system, Box 2.1. 

3.2.1 Adjustments to ADM equations and its effects on constraint propagations 

Generally, we can write the adjustment terms to (2.5) and (2.6) using (2.7) and (2.8) by the following 
combinations (using up to the first derivatives of constraints for simplicity): 



The adjusted ADM formulation [78]: Box 3.4 

Modify the evolution eqs of (jijjKij) by constraints H and Mi, 

adjustment terms of dtjij : +Pijn + Q^ijMk + P^ij{^kH) + q'^^jV^kMi), (3.7) 
adjustment terms of dtKij : +RijH + S'^ijMk + r''ij{Vk'H) + s''\j{VkMi), (3.8) 

where P, Q, R, S and p, q, r, s are multipliers. That is, 

daij = {2.h) + p^jH + Q^^jMk+p\j{Vk'n) + q^\,{VkMl), (3.9) 

dtKij = (2.Q) + Rijn + S^ijMk + r^ijiykn) + s^\,{VkMi). (3.10) 

According to this adjustment, the constraint propagation equations are also modified as 

dtH = (2.9) + Fr"(3.7) + Ht^^'dii?,.!) + Hi^'^''didj{2,.7) + i7r"(3.8), (3.11) 
dtMi = (2.10) + Mh"*"(3.7) + M2/'""aj(3.7) + M3i"*"(3.8) + Mu^'^''dj{2,.S). (3.12) 

with appropriate changes in indices. (See Appendix A. Hi, - ■ ■ , Mi, ■ ■ ■ are defined there.) 



We show two examples of adjustments here. We will list several others later in Table 2. 
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1. The standard ADM vs. original ADM 

The first comparison is to sliow tiie differences between tlie standard ADM [98] and the original 
ADM system [12] (see §2.0). In the notation of (3.7) and (3.8), the adjustment, 

Rij = KpO-'Jij, (and set the other multipliers zero) (3.13) 

will distinguish two, where Kp is a constant. Here Kp = corresponds to the standard ADM 
(no adjustment), and Kp = —1/4 to the original ADM (without any adjustment to the canonical 
formulation by ADM). As one can check by (3.11) and (3.12) adding Rij term keeps the constraint 
propagation in a first-order form. Prittelli [39] (see also [95]) pointed out that the hyperbolicity 
of constraint propagation equations is better in the standard ADM system. 

2. Detweiler type 

Dctwcilcr [33] found that with a particular combination, the evolution of the energy norm of the 
constraints, + Ai^, can be negative definite when we apply the maximal slicing condition, 
K = 0. His adjustment can be written in our notation in (3.7) and (3.8), as 



Pij = -KLa^lij, (3.14) 

Rij = KL(y''{Kij^{l/?,)K-i,j), (3.15) 

S\j = KLa'[?,{d^ia)5';) - {diaHj^% (3.16) 

s^\j = KLa%5^^)-{l/3Hjj''% (3.17) 



everything else is zero, where kl is a constant. Detweiler's adjustment, (3.14)-(3.17), does not 
put constraint propagation equation to a first order form, so we cannot discuss hyperbolicity 
or the characteristic speed of the constraints. We confirmed numerically, using perturbation on 
Minkowskii and Schwarzschild spacetime, that Detweiler's system provides better accuracy than 
the standard ADM, but only for small positive kl- 



3.2.2 Procedure to evaluate constraint amplification factors in spherically symmetric 
spacetime 

Before we compare between particular adjustment examples, we describe our procedure to evaluate 
CAFs in spherically symmetric spacetime. According to our motivation, the actual procedure to 
analyze the adjustments is to substitute the perturbed metric to the (adjusted) evolution equations first 
and to evaluate the according perturbative errors in the (adjusted) constraint propagation equations. 
However, for the simplicity, we apply the perturbation to the pair of constraints directly and analyze 
the effects of adjustments in its propagation equations. The latter, we think, presents the feature of 
constraint propagation more clearly for our purposes. 

The discussion becomes clear if we expand the constraint Cfj, := {H,Mi)^ using vector harmonics, 

C^ = J2 (^'"«^- + B^^bim + C'^Cim + D^'^dim) , (3.18) 
l,m 

where we choose the basis as 

aimie,ip) = iYi„„0,0,Of, (3.19) 
bim{9,^) = {0,Yi^,0,Of, (3.20) 
cim(.e,ip) = -j=^=={O,O,d0Yi^,d^Yi^f, (3.21) 

dim{e,^) = J—^ {0, 0, -J-d^Yim, sine deYimf, (3.22) 
Vi(i + 1) smt' 
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(a) 



(b) 



no adjustments (standard ADM) 



original ADH/I (k = - 1/4) 



Ol -0.5 

e 





Figure 5: Amplification factors (CAFs, eigenvalues of homogenized constraint propagation equations) 
are shown for the standard Schwarzschild coordinate, with (a) no adjustments, i.e., standard ADM, 
and (b) original ADM (kf = —1/4). [see eq. (3.13)]. The solid lines and the dotted lines with 
circles are real parts and imaginary parts, respectively. They are four lines each, but actually the two 
eigenvalues are zero for all cases. Plotting range is 2 < r < 20 using Schwarzschild radial coordinate. 
We set A; = 1,/ = 2, and m = 2 throughout the article. (Reprinted from [78], ©APS 2002) 



and the coefficients A , • • • , D are functions of {t, r). Here is the spherical harmonic function. 



Y, 



Im 



[OM = (-l)(-+H)/2 + 1) (; P,"^(cosg)e-^. 



The basis (3.19)-(3.22) are normalized so that they satisfy 



Jo Jo 



10 Jo 

where rj^^ is Minkowskii metric and the asterisk denotes the complex conjugate. Therefore 



A 



Im I „lm 



(a|-),a), dtA'^ = {a'C;^,dtC,), etc. 



(3.23) 



(3.24) 



(3.25) 



In order to analyze the radial dependences, we also express these evolution equations using the 
Fourier expansion on the radial coordinate. 



^Zm^^^Zm^^^gifcr ^^^^ 



(3.26) 



So that we can obtain the RHS of the evolution equations for (j4|™^(t), • • • , Z)|^(t))^ in a homogeneous 
form. 



3.2.3 Constraint amplification factors in Schwarzschild spacetime 

We present our CAF analysis in Schwarzschild black hole spacetime, which metric is 



2M, , 9 dr^ 

-(1 )dt^ + T- 

^ r ' 1 - 2M r 



+ r dO. ^ (the standard expression) 



(3.27) 



where M is the mass of a black hole. For numerical relativists, evolving a single black hole is the 
essential test problem, though it is a trivial at first sight. The standard expression, (3.27), has a 
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5 10 15 20 5 10 15 20 



Figure 6: Amplification factors of the standard Schwarzschild coordinate, with Detweiler type adjust- 
ments, (3.14)-(3.17). Multiphers used in the plot are (a) = +1/2, and (b) k,l = —1/2- Plotting 
details are the same as Fig.5. (Reprinted from [78], ©APS 2002) 

coordinate singularity at r = 2M, so that we need to move another coordinate for actual numerical 
time integrations. 

The ingoing Eddington-Finkelstein (iEF) coordinate has become popular in numerical relativity, in 
order to excise black hole singularity, since iEF penetrates the horizon without an irregular coordinate. 
The expression is, 

ds^ = -(1 - —)dtfp;p + —dtiEpdr + (1 + —)dr'^ + r^dO^ (the iEF expression) (3.28) 

which is given by tiEF = t + 2Mlog(r — 2M) and the radial coordinate is common to (3.27). 
We show CAFs in the following cases. More examples are available in [78]. 

1. in the standard Schvi^arzschild metric expression (3.27) 

• For the adjustment (3.13), CAFs are obtained as 

A' = (0,0, Va,-\/«), (3.29) 

2 4:Mk'^r'^{r-M) + 2M{2r-M) + l{l + l)r{r-2M) + ikr{2r^-Udr-2M'^) 
a - -k + ^ 

for the choice of ki = (the standard ADM), while they are 

. / /T /7n , M(2r - M) + irkM(2M - r) , , 

A* = 0,0,\/6,-\/6 , h=^ ^ '- 3.30 

for the choice of ki = —1/4 (the original ADM). 

These are plotted in Fig.5. The solid lines and dotted lines with circles are real parts and 
imaginary parts of CAFs, respectively. They are four lines each, but as we showed in (3.29), 
two of them are zero. The plotting range is 2 < r < 20 in Schwarzschild radial coordinate. 
The CAFs at r = 2 are zt-^/S/S and 0. The existence of this positive real CAF near the 
horizon is an important result. We show only the cases with I = 2 and k = \, because we 
judged that the plots of / = and other ks are qualitatively the same. 
The adjustment (3.13) with up = —1/4 returns the system back to the original ADM. 
CAFs are (3.30) and we plot them in Fig. 5(b). We can see that the imaginary parts are 
apparently different from those of the standard ADM [Fig. 5(a)]. This is the same feature as 
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Figure 7: CAFs in the iEF coordinate (3.28) on t = slice for the standard ADM formulation (i.e. no 
adjustments) [See Fig. 5(a) for the standard Schwarzschild coordinate.] and for Detweiler adjustments 
with kl = +1/2 [See Fig.6(a) for the standard Schwarzschild coordinate.]. The solid four lines and 
the dotted four lines with circles are real parts and imaginary parts, respectively. (Reprinted from 
[78], ©APS 2002) 

in the case of the flat background [95]. According to our conjecture, the non-zero imaginary 
values are better than zeros, so we expect that the standard ADM has a better evolution 
property than the original ADM system. Negative Kp always makes the asymptotical real 
values finite. 

• The Detweiler-type adjustment (3.14)-(3.17) makes the feature completely different. 
Fig. 6(a) and (b) are the cases of kl = ±1/2. A great improvement can be seen in the 
positive kl case where all real parts become negative in large r. Moreover all imaginary 
parts are apart from zero. These are the desired features according to our conjecture. 
Therefore we expect the Detweiler adjustment has good stability properties except near the 
black hole region. The CAF near the horizon has a positive real component. This is not 
contradictory with the Detweiler's original idea. His idea came from suppressing the total 
L2 norm of constraints on the spatial slice, while our plot indicates the existence of a local 
violation mode. The change of signature of kl can be understood just by changing the 
signature of CAFs, and this fact can also be seen to the other plot. In [78] we reported 
that a partial adjustment (apply only (3.14) or (3.17)) is also effective. 

2. in the iEF coordinates (3.28) 

• For the adjustment (3.13), we plotted CAFs in Figure 7. We Figure 7(a) is qualitatively 
different from Fig. 5(a). This is because the iEF expression is asymmetric to time, i.e. has 
non-zero extrinsic curvature. We notice that while some CAFs in iEF remain positive in 
large r region, that their nature changes due to the adjustments. 

• For the Detvi^eiler-type adjustment (3. 14)- (3. 17), CAFs are as in Figure 7(b). Inter- 
estingly, all plots indicate that all real parts of CAFs are negative, and imaginary parts are 
non-zero (again except near the black hole region). By arranging the multiplier parameter, 
there is a chance to get all negative real CAFs outside the black hole horizon. For example, 
all the real-part goes negative outside the black hole horizon \i hl > 3.1, while large kl 
may introduce another instability problem [94]. 

Such kinds of test can be done with other combinations. In Table 2, we listed our results for more 
examples. We defined the adjustment terms so that their positive multiplier parameter. At > 0, makes 
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No. 


adjustment 


1st? 


TRS 


Sch coord, 
real. 


imag. 


iEF coorc 
real. 


imag. 





— no adjustments 


yes 












P-1 


-Pit — KT Ot^^'i'i 


no 


no 


makes 2 Neg. 


not app. 


makes 2 Neg. 


not app. 


P-2 


Pi A — K r Q/yi -i 


no 


no 


makes 2 Neg. 


not app. 


maJjes 2 Neg. 


not app. 


P-3 




no 


no 


slightly enl. Neg. 


not app. 


slightly enl. Neg. 


not app. 


P-4 


P — K'V 

ij fij 


no 


no 


makes 2 Neg. 


not app. 


makes 2 Neg. 


not app. 


P-5 


^P^ j r 


no 


no 


red. Pos./enl.Neg. 


not app. 


red. Pos. /enl. Neg. 


not app. 


Q-1 


LJ ■ 1 ''J 


no 


no 


N/A 


N/A 


K ~ 1.35 min. vals. 


not app. 


Q-2 


2,j T"P ^ 


no 


yes 


red. abs vals. 


not app. 


red. abs vals. 


not app. 


Q-3 


Q ij Q ij '^Tij 


no 


yes 


red. abs vals. 


not app. 


cnl.Neg. 


enl. vals. 


Q-4 


ij 7*7" — f^^rf 


no 


yes 


red. abs vals. 


not app. 


red. abs vals. 


not app. 


R-1 


Rij KpOfyij 


yes 


yes 


kf = —1/4 min. abs vals. 


Kf = —1/4 min. vals. 


R-2 


R'ij R-pv — f^jj,^ 


yes 


no 


not app. 


not app. 


red.Pos./enl.Neg. 


enl. vals. 


R-3 


Rij Rrr — K'Yrr 


yes 


no 


enl. vals. 


not app. 


red . Pos . /enl . Neg. 


enl. vals. 


S-1 


S^ij S\j = (3.16) 


yes 


no 


not app. 


not app. 


not app. 


not app. 


S-2 




yes 


no 


makes 2 Neg. 


not app. 


makes 2 Neg. 


not app. 


p-1 


p'^ij p'"ij = -Kajij 


no 


no 


red. Pos. 


red. vals. 


red. Pos. 


enl. vals. 


p-2 


P^ij P^rr = Ka 


no 


no 


red. Pos. 


red. vals. 


red . Pos/enl . Neg. 


enl. vals. 


p-3 


P ij P rr — KOi'^rr 


no 


no 


makes 2 Neg. 


enl. vals. 


red.Pos.vals. 


red.vals. 


q-1 


q^Uj q^'^'a = «Q;7i3 


no 


no 


K = 1/2 min. vals. 


red.vals. 


not app. 


enl. vals. 


q-2 


Q. ij rr — I^Oi'~^rr 


no 


yes 


red. abs vals. 


not app. 


not app. 


not app. 


r-1 


T ij T ij = KOC^ij 


no 


yes 


not app. 


not app. 


not app. 


enl. vals. 


r-2 


^ ij ^ rr — tZOt 


no 


yes 


red. abs vals. 


enl. vals. 


red. abs vals. 


enl. vals. 


r-3 


V \j V rr — KOL'yrr 


no 


yes 


red. abs vals. 


enl. vals. 


red. abs vals. 


enl. vals. 


s-1 


S ij S ij = (3.17) 


no 


no 


makes 4 Neg. 


not app. 


makes 4 Neg. 


not app. 


s-2 


^ ij ^ ij — f^f^'Yij 


no 


no 


makes 2 Neg. 


red. vals. 


makes 2 Neg. 


red. vals. 


s-3 


S 'ij S rv — KCX'Yrv 


no 


no 


makes 2 Neg. 


red. vals. 


makes 2 Neg. 


red. vals. 



Table 2: List of ADM adjustments we tested in the Schwarzschild spacetime. The column of adjust- 
ments are nonzero multipliers in terms of (3.7) and (3.8). The column 'TRS' indicates whether each 
adjusting term satisfies the time reversal symmetry or not on the standard Schwarzschild coordinate. 
('No' means a candidate that makes asymmetric CAFs.) The column '1st?' indicates whether each 
adjusting term breaks the first-order feature of the standard constraint propagation equation, (2.9) 
and (2.10). ('Yes' keeps the system first-order, 'No' may break hyperbolicity of constraint propagation 
depending a choice of k.) The effects to CAFs (when k > 0) are commented for both coordinates and 
both real/imaginary parts, respectively. The 'N/A' means that there is no effect due to the coordinate 
properties; 'not app.' (not apparent) means the adjustment does not change the CAFs effectively 
according to our conjecture; 'enl. /red. /min.' means enlarge/reduce/minimize, and 'Pos. /Neg.' means 
positive/negative, respectively. These judgements are made at the r ^ O(IOM) region on their t = 
slice. See more detail in [78]. 

the system better in stability according to our conjecture. (Here better means in accordance with our 
conjecture in Box 3.2 and 3.3). The table includes the above results and is intended to extract the 
contributions of each term in (3.7) and (3.8). The effects of adjustments (of each k > case) to CAFs 
are commented upon for each coordinate system and for real/imaginary parts of CAFs, respectively. 
These judgements arc made at the r ~ O{10M) region on their t = slice. Among them, No. R-2 
in Table 2 explains why a particular adjustment by PennState group [51] gives better stability than 
before. 

3.2.4 Remarks 

Numerical demonstrations The above analyses are only predictions, and supporting numerical 
demonstrations are necessary for the next steps. Systematic numerical comparisons are progressing, 

and we show two sample plots here. 

Figure 8 (a) is a test numerical evolution of Detweiler-type adjustment on the Minkowskii back- 
ground. We see the adjusted version gives convergence on to the constraint surface by arranging the 



28 



ADM and its adjusted versions 




300 



time time 



Figure 8: Comparisons of numerical evolution between adjusted ADM systems, (a) Demonstration of 
the Detweiler's modified ADM system on Minkowskii background spacetime, 1-dimensional simulation. 
The L2 norm of the constraints fi^^^'^ and A4^^^^ is plotted in the function of time. Artificial 
error was added at t = 0.25. L is the parameter used in (3.14)-(3.16). We see the evolution is 
asymptotically constrained for small k > 0. (Reprinted from [95], ©APS 2001) (b) L2 norm of the 
Hamiltonian constraint 7^^^^^ of evolution using ADM/adjusted ADM formulations for the case of 
Teukolsky wave, 3-dimensional simulation. Cactus-based original ADM module (CactusGR) was used. 
(Shinkai-Yoneda, in preparation). 



magnitude of the adjusting parameter, k. 

Figure 8 (b) is an example from the project of "Comparison of formulations of the Einstein equa- 
tions for numerical relativity" [59]. (We will describe the project in §4). The plot was obtained by a 
3-dimensional numerical evolution of weak gravitational wave, the so-called Teukolsky wave [85]. The 
lines are of the original/standard ADM evolution equations, Detweiler-type adjustment, and a part of 
Detweiler-type adjustment (actually used only (3.14)). For a particular choice of n, we observe again 
the L2 norm of constraint (violation of constraints) is reduced than the standard ADM case, and can 
evolve longer than that. 

Notion of Time Reversal Symmetry During the comparisons of adjustments, we found that it 
is necessary to create time asymmetric structure of evolution equations in order to force the evolution 
on to the constraint surface. There are infinite ways of adjusting equations, but we found that if we 
follow the guideline Box 3.5, then such an adjustment will give us time asymmetric evolution. 



Trick to obtain asymptotically constrained system: Box 3.5 

= Break the time reversal symmetry (TRS) of the evolution equation. 

1. Evaluate the parity of the evolution equation. 

By reversing the time {dt — > —dt), there are variables which change their signatures (parity 
(-)) [e.g. Kij, dt'jij, Mi, ■■■], while not (parity (+)) [e.g. gij,dtKij,7i,- ■ ■]■ 

2. Add adjustments which have different parity of that equation. 

For example, for the parity (— ) equation dt^ij, add a parity (+) adjustment kH. 
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One of our criteria, the negative real CAFs, requires breaking the time-symmetric features of the 
original evolution equations. Such CAFs are obtained by adjusting the terms which break the TRS of 
the evolution equations, and this is available even at the standard ADM system. The TRS features (on 
the standard Schwarzschild metric, which is time-symmetric background) are also denoted in Table 2. 
The criteria of Box 3.5 will be applied also to the BSSN equation (see §3.3.3). 

Differences with Detweiler's requirement We comment on the differences between Detweiler's 
criteria [33] and ours. Detweiler calculated the L2 norm of the constraints, Cp, over the 3-hypersurface 
and imposed the negative definiteness of its evolution, 

Detweiler's criteria ^ dt j CpC dV V non zero Cp. (3.31) 

where CpC^ =: G'^'^CpCa, and Gpa = diag[l,^ij] for the pair of Cp = {H,Mi). 

Assuming the constraint propagation to be dfCp = Ap'^C^ in the Fourier components, the time 
derivative of the L2 norm can be written as 

dtiCpC") = {AP" + A^'P + dtGP'')Cpha. (3.32) 

Together with the fact that the L2 norm is preserved by Fourier transform, we can say, for the case 
of static background metric, 

Detweiler's criteria ^ eigenvalues of (^ + J^) are all negative Vfc. (3.33) 
Our criteria eigenvalues of A are all negative V/c. (3.34) 

Therefore for the case of static background, Detweiler's criterion is stronger than ours. For example, 
the matrix 

A= where o is constant, (3.35) 

for the evolution system (C\^C2) satisfies our criterion but not Detweiler's when \a\ > \pl. This 
matrix however gives asymptotical decay for ((71,(72). Therefore we may say that Detweiler requires 
the monotonic decay of the constraints, while we assume only asymptotical decay. 

We remark that Detweiler's truncations on higher order terms in (7-norm corresponds to our 
perturbational analysis; both are based on the idea that the deviations from constraint surface (the 
errors expressed non-zero constraint value) are initially small. 

Ranges of effective k We do not discuss the ranges of the effective multiplier parameter, k, since 

the range depends on the characteristic speeds of the models and numerical integration schemes as we 
observed in [94]. At this moment, we can only estimate the maximum value for adjusting parameter 
KS from von Neumann's stability analysis, while we do not have a background theoretical explanation 
to predict the optimized value of k for numerical evolution. We will comment this point later again 
in §4. 

3.3 Adjusted BSSN formulations 

Next we apply the idea of adjusted system to the BSSN system. Box 2.3. 
3.3.1 Constraint propagation analysis of the BSSN equations 

The BSSN has 5 constraint equations, (W, Mi, Q\ A, 5), see (2.26)-(2.30). We begin identifying where 
the BSSN evolution equations were adjusted already in its standard notation, (2.21)-(2.25). 
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Taking careful account of these constraints, (2.26) and (2.27) can be expressed directly as 

U = e-^'^R-8e-'^'PDWjip-8e-^'P{&^p){Dj^p) + {2/3)K^ -AijA'^ -{2/3)AK, (3.36) 
Mi = 6A^i{Djcp)-2A{Diip)-{2/3){DiK)+j>'^{DjAki). (3.37) 

By a straightforward calculation, we get: 

dfcp = d^cp + {l/6)aA-{l/12)r\djS)P^, (3.38) 

dhij = dfyi,-{2/3)aji,A + {l/3)r\dkS)P%j, (3.39) 
d^K = dfK-{2/3)aKA-an + ae-*'P{Djg^), (3.40) 
d^Aij = d^Aj + {{l/3)a%K -i2/3)aAi,)A 

+{il/2)ae-^^i^k7^J) - {l/6)ae'^^^^,r\dkS))g'' 

+«e-4'^7fc(.(5,)^') - (l/3)ae-^'^7,,(a,.g^), (3.41) 
r = d^r + ( - (2/3)(9,a)f ^ - i2/3)a{d,r) - {l/3)arr\d,S) + iaf^d,ip))A 

~{2/3)af\d,A) + 2aY^Mj - {l/2){dk/3^^)^''^rHdjS) + {l/6){dj (3'') f^r\dkS) 
+{l/3){dkP'')f^rHd,S) + i5/6)p'^r'r{dkS)idjS) + {l/2)(5^r\dkrmS) 
+[^/3)(3^r\dnn{dkS). (3.42) 

where denotes the part of no replacements, i.e. the terms only use the standard ADM evolution 
equations in its time derivatives. 

Prom (3.38)-(3.42), we understand that all BSSN evolution equations are adjusted using constraints. 
However, from the viewpoint of time reversal symmetry (Box 3.5), all the above adjustments in (3.38)- 
(3.42) unfortunately keep the time reversal symmetry. Therefore we can not expect direct decays of 
constraint violation in the present form. 

The set of the constraint propagation equations, dtCH, Aii,Q^ , A,S)'^ , turns to be not a first-order 
hyperbolic and includes many non-linear terms, (see the Appendix in [96]). In order to understand 
the fundamental structure, we show an analysis on the flat spacetime background. 

For the flat background metric g^i, = 1)^1,, the first order perturbation equations of (3.38)-(3.42) 
can be written as 

dt% = -{l/6pK+{l/6){K^-lpA, (3.43) 
dt^%j = -i^^ij - {2/3)iKy - l)5ij'^^U, (3.44) 
dt^^^K = -{djdl% + {KKi-l)dj'^^^^ -iKK2-lf^, (3.45) 
= ('\Rff^''f^ - ^'\DiDjaf^ + {KAi - l)5u^i{d^j'^^) - {1/3){KA2 - l)%(5,(i^t^46) 
a/i)f^ = -(4/3)(a/ik) - (2/3) (Kfi - \){dPA) + 2(/Cf2 - IpMi, (3.47) 

where we introduced parameters ks, all k = reproduce no adjustment case from the standard ADM 
equations, and all k = 1 correspond to the BSSN equations. We express them as 

i^adj '■= {K'<p,Kj,kki,kk2,kai,ka2, 1^1,1^2)- (3.48) 
Constraint propagation equations at the first order in the flat spacetime, then, become: 

= (k^ - (2/3)Kp^ - (4/3)k^ + 2) djd/^U + 2(Kp2 - l){dj'^^^Mj), (3.49) 
dt^^^Mi = i-{2/3)KKi + {1/2)kai - (1/3)ka2 + (1/2)) didl^^^ 

+{l/2)f,Aidjdl''^' + {{2/3)kk2 - (1/2)) (3.50) 

dl'^' = 2Kf «M + (-(2/3)Kp,-(l/3)Av;^)(a/^U), (3.51) 

dt^^h = -2k4^U, (3.52) 

dl'U = {KAi-KA2)idl'^n- (3-53) 
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3.3.2 The origin of the advantages of the BSSN equations 

We next discuss CAFs of (3.49)-(3.53). Hereafter we let k"^ = k1 + ky + k1 for Fourier wave numbers. 

1. The no-adjustment case, Kadj =(all zeros). This is the starting point of the discussion. In this 
case, 

CAFs = (0(x7),±y-A:2), 

i.e., (0 (x7), zbpure imaginary (1 pair)). In the standard ADM formulation, which uses {'yij,Kij), 
CAFs are (0, 0, itPure Imaginary) [95]. Therefore if we do not apply adjustments in BSSN 
equations the constraint propagation structure is quite similar to that of the standard ADM. 

2. For the BSSN equations, Kadj =(all Is), 

CAFs = (0(x3),±V'^(3 pairs)), 

i.e., (0 (x3), ibPure Imaginary (3 pairs)). The number of pure imaginary CAFs is increased over 
that of No.l, and we conclude this is the advantage of adjustments used in BSSN equations. 

3. No 5-adjustment case. All the numerical experiments so far apply the scaling condition S for 
the conformal factor cp. The 5-originated terms appear many places in BSSN equations (2.21)- 
(2.25), so that we guess non-zero iS is a kind of source of the constraint violation. However, since 
all (S-originated terms do not appear in the flat spacetime background analysis, [no adjusted 
terms in (3.43)-(3.47)], our case can not say any effects due to S-constraint. 

4. No .4-adjustment case. The trace (or traceout) condition for the variables is also considered 
necessary (e.g. [3]). This can be checked with Kadj = {i^, f^, 1; Ij 1; 1) i^, 1), and we get 

CAFs = (0(x3),±v/^(3 pairs)), 
independent of k. Therefore the effect of ^-adjustment is not apparent from this analysis. 

5. No t?*- adjustment case. The introduction of F* is the key in the BSSN system. If we do not 
apply adjustments by Q\ {nadj = (1; 1)0, 1,0, 0, 1, 1)) then we get 

CAFs = (0(x7),±V^), 

which is the same with No.l. That is, adjustments due to Q'^ terms are effective to make a 
progress from ADM. 

6. No Alj-adjustment case. This can be checked with Kadj = (Ij 1) Ij 1) Ij 1) Ij z^)) and we get 

CAFs = (0,±V-AcA;2 (2 pairs), 

±-/^F^^^lT4^cT^^^^4^V6, ±y^P^^lT4^^^^^4^^0/6). 

If K = 0, then (0(x7), i^/k'^^/S), which is (0(x7), ibrcal value). Interestingly, these real values 
indicate the existence of the error growing mode together with the decaying mode. Alcubierre 
et al. [5] found that the adjustment due to the momentum constraint is crucial for obtaining 
stability. We think that they picked up this error growing mode. Fortunately at the BSSN limit 
{k = 1), this error growing mode disappears and turns into a propagation mode. 

7. No 7^-adjustment case. The set Kadj = (1) Ij 1) 1) 1) Ij 1) gives 

CAFs = (0(x3),±\/^(3 pairs)), 
independently to k. Therefore the effect of 7i-adjustment is not apparent from this analysis. 
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No. 




Constraints (number of components) 




diag? 


CAFs 






n (1) 


Mi (3) 


(3) 


-4(1) 


5 (1) 




in Minkowskii background 


0. 


standard ADM 


use 


use 


- 


- 




yes 


(0, 0, 


1. 


BSSN no adjustment 


use 


use 


use 


use 


use 


yes 


(0,0,0,0,0,0,0,0,0) 


2. 


the BSSN 


use+adj 


use+adj 


use+adj 


use+adj 


use+adj 


no 


(0, 0, 0, 0, 0, 0, 0, 0, 0) 


3. 


no 5 adjustment 


use+adj 


use+adj 


use+adj 


use+adj 


use 


no 


no diflference in flat background 


4. 


no A adjustment 


use+adj 


use+adj 


use+adj 


use 


use+adj 


no 


(0,0, 0,9;, Of, Of, Of, Of, Of) 


5. 


no 6' adjustment 


use+adj 


use+adj 


use 


use+adj 


use+adj 


no 


(0,0,0, 0, 0, 0,0,0, 0) 


6. 


no Mi adjustment 


use+adj 


use 


use+adj 


use+adj 


use+adj 


no 


(0,0,0,0,0, 0,0, SR, SR) 


7. 


no Ti adjustment 


use 


use+adj 


use+adj 


use+adj 


use+adj 


no 


(0,0, 0,0,0,0, 0, 0, 0) 



Table 3: Contributions of adjustments terms and effects of introductions of new constraints in the 
BSSN system. The center column indicates whether each constraints are taken as a component of 
constraints in each constraint propagation analysis ('use'), and whether each adjustments are on 
('adj'). The column 'diag?' indicates diagonalizability of the constraint propagation matrix. The 
right column shows amplification factors, where and means pure imaginary and real eigenvalue, 
respectively. 

These tests are on the effects of adjustments that are already in the BSSN equations. We will consider 
whether much better adjustments are possible in the next section. 

We list the above results in Table 3. The most characteristic points of the above arc No. 5 and 
No. 6 that denote the contribution of the momentum constraint adjustment and the importance of 
the new variable r\ It is quite interesting that the unadjusted BSSN equations (case 2) does not 
have apparent advantages from the ADM system. As we showed in the case 5 and 6, if we missed a 
particular adjustment, then the expected stability behaviour occationally gets worse than the starting 
ADM system. Therefore we conclude that the better stability of the BSSN formulation is obtained by 
their adjustments in the equations, and the combination of the adjustments is in a good balance. 

3.3.3 Proposals of the modified BSSN equations 

We next consider the possibility whether we can obtain a system which has much better properties; 
whether more pure imaginary CAFs or negative real CAFs. 

Heuristic examples 

(A) A system which has 8 pure imaginary CAFs: 

One direction is to seek a set of equations which make fewer zero CAFs than the standard BSSN 
case. Using the same set of adjustments in (3.43)-(3.47), CAFs are written as 

CAFs = ^0, ±^ —k'^KAiKf,2 (2 pairs), ibcomplicated expression, itcomplicated expression 

The terms in the first line certainly give four pure imaginary CAFs (two positive and negative 
real pairs) if KyiiKpg > (< 0). Keeping this in mind, by choosing Kadj = (1, 1, Ij 1, 1, 1"^, 1, 1), we 
find 

CAFs = (0,±y/^ (2 pairs), ±^-fe2(2 + k + \k- 4|)/6, ±^J-k^{2 + k - \k - 4|)/6, ) . 

Therefore the adjustment Kadj = (1, 1; 1, 1; 1,4, 1, 1) gives CAFs = (O, (4 pairs)) , which 

is one step advanced from BSSN according our guidelines. 

We note that such a system can be obtained in many ways, e.g. Kadj = (0, 0, 1, 0, 2, 1, 0, 1/2) also 
gives four pairs of pure imaginary CAFs. 
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(B) A system which has negative real CAF: 

One criterion to obtain a decaying constraint mode (i.e. an asymptotically constrained system) 
is to adjust an evolution equation as it breaks time reversal symmetry (Box 3.5). For example, 
we consider an additional adjustment to the BSSN equation as 

dtlij = dtlij + i^SDOLiijH, (3.54) 

which is a similar adjustment of the simplified Detweiler-type [33]. The constaint amplification 
factors become 

CAFs = (0 (x2), ±7^(3 pairs), (3/2)A:2Ksi5), 
in which the last one becomes negative real if ksd < 0. 

(C) Combination of above (A) and (B) 
Naturally we next consider both adjustments: 

dtlij = df^ij + KsDaiijTi (3.55) 
dtA, = dfAij- K^ae-^^^ijdkQ^ (3.56) 

where the second one produces the 8 pure imaginary CAFs. We then obtain 

CAFs = (0,±\/^(3 pairs), (3/4)A;^K5d ± \/k'^{-K8 + {9 / 16) k"^ ksd)) 

which reproduces case (A) when k^x) = 0, = 1, and case (B) when Kg = 0. These CAFs can 
become (0, pure imaginary (3 pairs), complex numbers with a negative real part (1 pair)), with 
an appropriate combination of Kg and ksd- 



Possible adjustments In order to break time reversal symmetry of the evolution equations (Box 
3.5), the possible simple adjustments are (1) to add 7i, S or terms to the equations of dt4>, dtlij, 
or dtT^, and/or (2) to add Mi or A terms to dtK or dtAij. We write them generally, including the 
above proposal (B), as 

dt(t) = df(j) + K^H(M + K^abkg^, (3.57) 
dtlij = df^ij + KSD a^ijU + KjQi a^ijbkQ^ + K^g2 a%{ibj-^g^ 

+ K^si otiijS + /«^52 aDiDjS, (3.58) 

dtK = dfK + KKMa¥\DjMk), (3.59) 
dtAij = dfAij + KAMI oiiij {b^Mk) + KAM2 a{b^iMj) ) 

+ K-AAi O-iijA + KAA2 abibjA, (3.60) 

= dfr + Kf^a&n + Kfg^ag' + Kfgr^abWjg^ + Kfg^a&bjg^, (3.61) 

where ks are possible multipliers (all k = reduce the system the standard BSSN system). 

Wc show the effects of each terms in Tabic 4. The CAFs in the table are on the fiat space 
background. We see several terms make negative CAFs, which might improve the stability than the 
previous system. For the readers convenience, we list up several best candidates here. 

(D) A system which has 7 negative CAFs 

Simply adding D(^iMj) term to dtA^j equation, say 

dtAij = 5f ^'^^^Ay + KAM2a{b^iMj)) (3.62) 
with KAM2 > 0, CAFs on the fiat background are 7 negative real CAFs. 
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adjustment 


CAFs 


diag? 


cfTect of the adjustment 


dt4> 




(0,0, ±\/ 


-fc2(*3),8K^Hfe^) 


no 


i^4>'H. < makes 1 Ncg. 




dt<l> 




(0,0, ±^ 


— fc2(*2), long expressions) 


yes 


>^<l>g < makes 2 Neg. 1 Pos. 




dtlij 


KSD CXlijTi 


(0,0, ±\/ 


^(*3), (3/2)Ksi3A;2) 


yes 


KSD < makes 1 Neg. 


Case (B) 


dtlij 


^751 a^ijDkQ^ 


(0,0, ±v^ 


— fe^(*2), long expressions) 


yes 


K^gi > makes 1 Neg. 




dtlij 


^752 aik{iDj)g'' 


(0,0, long 


expressions) 


yes 


«^7e2 < makes 6 Neg. 1 Pos. 


Case (El) 


dtlij 


K7S1 a^ijS 


(0,0, ±/ 


=^(*3),3k^si) 


no 


KXfSi < makes 1 Neg. 




dtlij 


i^^S2 OiDiDjS 


(0,0, ±^ 


-fc2(*3),-K-52fc2) 


no 


"752 > makes 1 Neg. 




dtK 


liKMOil^^iDjMk) 


(0, 0,0,± 


V — fc2(*2), long expressions) 


no 


i^KM < makes 2 Neg. 




dtAjj 


KAMI 0'lij{D''Mk) 


(0,0, iV" 


-fc2(*3), -KAMlfc^) 


yes 


i^AMi > makes 1 Neg. 






KAM2 0:{D(iMj)) 


(0,0, long 


expressions) 


yes 


i^AM2 > makes 7 Neg 


Case (D) 


dtAij 


I^AAl O^lijA 


(0,0, ±\/ 


-A;2(*3),3ka^i) 


yes 


i^AAi < makes 1 Neg. 




dtAij 


K-AA2 aDiDjA 


(0,0, ±^ 


-k^ -K.AA2k^) 


yes 


'^AA2 > makes 1 Neg. 




dtf' 


aD'-'H 


(0,0, ±/ 


-fe2(*3), -KA^2fc^) 


no 


'^r'H ^ ^ makes 1 Neg. 




dtf* 




(0,0, long 


expressions) 


yes 


'^rgi ^ ^ makes 6 Neg. 1 Pos. 


Case (E2) 


dtV* 




(0,0, long expressions) 


yes 


Kpg2 > makes 2 Neg. 1 Pos. 




9tfi 




(0,0, long expressions) 


yes 


Kpug > makes 2 Neg. 1 Pos. 





Table 4: Possible adjustements which make a real-part constraint amplification factors negative. The 
column of adjustments are nonzero multipliers in terms of (3.57)-(3.61), which all violate time reversal 
symmetry of the equation. The column 'diag?' indicates diagonalizability of the constraint propagation 
matrix. Neg./Pos. means negative/positive respectively. 



(E) A system which has 6 negative and 1 positive CAFs 

The below two adjustments will make 6 negative real CAFs, while they also produce one positive 
real CAF (a constraint violating mode). The effectiveness is not clear at this moment, but we 
think they are worth to be tested in numerical experiments. 

(^1) dt^ij = d^^^^% + K^g2a^k(iDj)Q^, with K;yg2<0. (3.63) 

{E2) dtV = df^^^P + Kpg2 aDWjg\ with Kfg^ < 0. (3.64) 

3.3.4 Remarks 

We have studied step- by-step where the replacements in the equations affect and/or newly added 
constraints work, by checking whether the error of constraints (if it exists) will decay or propagate 
away. Alcubierre et al [5] pointed out the importance of the replacement (adjustment) of terms in 
the evolution equation due to the momentum constraint, and our analysis clearly explain why they 
concluded this is the key. Not only this adjustment, we found, but also other adjustments and other 
introductions of new constraints also contribute to making the evolution system more stable. We 
found that if we missed a particular adjustment, then the expected stability behaviour occationally 
gets worse than the ADM system. We further propose other adjustments of the set of equations which 
may have better features for numerical treatments. 

The discussion was only in the fiat background spacetime, and further analysis is in progress. 
However, we rather believe that the general fundamental aspects of constraint propagation analysis 
are already revealed here. This is because, for the ADM and its adjusted formulation cases, we found 
that the better formulations in the flat background are also better in the Schwarzschild spacetime, while 
there are differences on the effective adjusting multipliers or the effective coordinate ranges [78, 95]. 
Actually, recently Yo, Baumgarte and Shapiro [90] reported their simulations of stationary rotating 
black hole, and mentioned that the above proposal (B) was contributed to maintain their evolution of 
Kerr black hole (J/M up to 0.9M) for long time (t ~ 6000M). Their results also indicates that the 
evolved solution is closed to the exact one, that is, the constrained surface. 
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4 Outlook 



4.1 What we have achieved 

Let us summarize our story first. The beginnings of the study are: 

• General relativistic numerical simulations are quite important in astrophysical studies, but we do 
not have a definite recipe to integrate the Einstein equations in a long-term stable and accurate 
manner. 

• The most standard approach is to decompose the space-time into 3-space and time (3 + 1 
decomposition), to solve the constraints for obtaining the initial data, and to evolve the space- 
time applying a free evolution scheme. 

• Over a period of decades, the community has been observing the violation of constraints in their 
simulations, and their consensus is that we have to find a better formulation of the Einstein 
evolution equations. 

We reviewed recent efforts on this problem by categorizing them into 

(0) The standard ADM formulation (§2.0, Box 2.1), 

(1) The modified ADM (so-called BSSN) formulation (§2.1, Box 2.3), 

(2) Hyperbolic formulations (§2.2, Box 2.5), and 

(3) Asymptotically constrained formulations (§2.3). 

Among them, the approach (2) is perhaps best justified on mathematical grounds. However, as we 
critically reviewed in §2.2.3, the practical advantages may not be available unless we kill the lower-order 
terms in its hyperbolized equation, as in KST's formulation (Box 2.6). 

We therefore proceeded in the direction (3). Our approach, which we term adjusted system, is to 
construct a system that has its constraint surface as an attractor. Our unified view is to understand 
the evolution system by evaluating its constraint propagation. Especially we proposed to analyze the 
constraint amplification factors (Box 3.1) which are the eigenvalues of the homogenized constraint 
propagation equations. We analyzed the system based on our conjecture (Box 3.2/3.3) whether the 
constraint amplification factors suggest the constraint to decay /propagate or not. We concluded that 

• The constraint propagation features become different by simply adding constraint terms to the 
original evolution equations (we call this the adjustment of the evolution equations) . 

• There is a constraint-violating mode in the standard ADM evolution system when we apply it 
to a single non-rotating black hole space-time, and its growth rate is larger near the black-hole 
horizon. 

• Such a constraint-violating mode can be killed if we adjust the evolution equations with a 
particular modification using constraint terms (Box 2.7). An effective guideline is to adjust 
terms as they break the time-reversal symmetry of the equations (Box 3.5). 

• Our expectations are borne out in simple numerical experiments using the Maxwell, Ashtekar, 
and ADM systems. However the modifications are not yet perfect to prevent non-linear growth 
of the constraint violation. 

• We understand why the BSSN formulation works better than the ADM one in the limited 
case (perturbative analysis in the flat background), and further we proposed modified evolution 
equations along the lines of our previous procedure. 

The common key to the problem is how to adjust the evolution equations with constraints. Any 
adjusted systems are mathematically equivalent if the constraints are completely satisfied, but this is 
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not the case for numerical simulations. Replacing terms with constraints is one of the normal steps 
when people hyperbolize equations. Our approach is to employ the evaluation process of constraint 
amplification factors for an alternative guideline to hyperbolization of the system. 

4.2 Next steps 

Here are some directions for future researches in this formulation problem of the Einstein equations. 

Generalize the adjusted systems: We have tried to unify all the efforts of the community applying 
the idea of "adjusted systems" . Our newly modified equations are working as desired up to the current 
numerical tests, but we see also that the effect is not yet perfect to remove non-linear error growth 
that terminates numerical simulations. Wc suspect that this is due to the current eigenvalue analysis 
based on a perturbative analysis with fixing adjusting multiplier, n. One remedy is to generalize the 
determination process of k, say, dynamically and automatically under a suitable principle. We are 
now working on a method to control the violation of each constraint independently, or an additional 
supporting mathematical criteria to realize more robust stabilizations. 

More on hyperbolic formulations: We have already pointed out several directions on hyperbolic 
efforts in §2.2.3. We here only point out the links to the initial-boundary value problem (IB VP). In 
order to avoid unphysical incoming information from the boundary of computation, a proper treatment 
of the boundary is quite important problem in numerical simulations. If we can treat boundary 
condition as a part of mathematical framework, that would be a great step to promote researches. 
The IBVP approach to the Einstein equations is quite new and has only been studied in a restricted 
cases. The current proposals are based on a particular symmetric hyperbolic formulation or under a 
certain assumption to the symmetry of space-time. More general treatments on IBVP are expected. 

Alternative new ideas?: We have to be open to develop an alternative approach to this formulation 
problem. If our goal is to obtain a stable system on its constrained surface, then there may be at least 
three fundamental approaches: (1) construct an improved evolution system (as we discussed most), 
(2) develop a "maintenance method" for the system, or (3) redefine or classify the stability of the 
system. An idea of automatic adjustment of multipliers in adjusted system is a sort of "maintenance". 
We guess there might be a key in existing theories such as control theories, optimization methods 
(convex functional theories), mathematical programming methods, or others. 

Numerical comparisons of formulations: Fortunately, an effort toward systematic numerical 
comparisons of different formulations of the Einstein equations has been organized and started recently. 
The workshop "Comparisons of Formulations of Einstein's equations for Numerical Relativity" was 
held at Mexico City in May 2002 ■^ and more than 20 people attended this 2-wcck workshop. This 
project intends to provide a format for comparisons in a common numerical framework. That is, 
comparing different formulations using the same initial data, the same resolution, the same integration 
scheme, the same boundary treatment, and the same output. Comparisons are now in progress for 
vacuum and regular space-time, and we hope to extend them to black-hole space-time next. The first 
reports are in preparation [59]. 

4.3 Final remcirks 

If we say the final goal of this project is to find a robust algorithm to obtain long-term accurate 
and stable time-evolution method, then the recipe should be a combination of (a) formulations of the 

^Follow-up information is available at http://www.nuclecu. imam. mx/"'gravit/main_m. html. 
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evolution equations, (b) choice of gauge conditions, (c) treatment of boundary conditions, and (d) 
numerical integration methods. We are in the stages of solving this mixed puzzle. The ideal almighty 
algorithm may not exit, but we believe our accumulating experience will make the ones we do have 
more robust and automatic. 

We have written this review from the viewpoint that the general relativity is a constrained dynam- 
ical system. This is not only a proper problem in general relativity, but also in many physical systems 
such as electrodynamics, magnetohydrodynamics, molecular dynamics, mechanical dynamics, and so 
on. Therefore sharing the thoughts between different field will definitely accelerate the progress. 

When we discuss asymptotically constrained manifolds, wc implicitly assume that the dynamics 
could be expressed on a wider manifold. Recently we found that such a proposal is quite similar 
with techniques in e.g. molecular dynamics. For example, people assume an extended environment 
when simulating molecular dynamics under constant pressure (with a potential piston) or a constant 
temperature (with a potential thermostat) (see e.g. [65]). Wc have also noticed that a dynamical 
adjusting method of Lagrange multipliers has been developed in multi-body mechanical dynamics (see 
e.g. [62]). We are now trying to apply these ideas back into numerical relativity. 

In such a way, communication and interaction between different fields is encouraged. Cooperation 
between numerical and mathematical scientists is necessary. By interchanging ideas, we hope we will 
reach our goal in next few years, and obtain interesting physical results and predictions. It is our 
personal view that exciting revolutions in numerical relativity are coming soon. 
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A General expressions of ADM constraint propagation equations 

For the reader's convenience, we express here the constraint propagation equations generahy, consid- 
ering the adjustments to the evolution equations. 

A.l The standard ADM equations and constraint propagations 

We start by analyzing the standard ADM system, that is, with evolution equations (2.5) and (2.6) 
and constraint equations (2.7) and (2.8). 

The constraint propagation equations, which are the time evolution equations of the Hamiltonian 
constraint (2.7) and the momentum constraints (2.8). 

Expression using H and Mi The constraint propagation equations can be written as [these are 
the same with (2.9) and (2.10) ] 

dtU = {djH) + 2aKn - 2a-f'^ (diMj) 

+a(az7™fc)(27™'7'^' - 7™S'^)-M, - iY'{^Ja)M^, (A.l) 
dtMi = -{l/2)a{din)-{dia)rL + (3^{d,Mi) 

+aKMi - P^-f^\daik)Mj + {dA)^''^Mj. (A.2) 

This is a suitable form to discuss hyperbolicity of the system. The simplest derivation of (A.l) and 
(A.2) is by using the Bianchi identity, which can be seen in Prittelli [39]. 
A shorter expression is available, e.g. 

dtU = /3^din + 2aKn-2a-f-^/^di{^M^)-A{dia)M^ 

= P^ViH + 2aKH-2a{ViM^)-A{Via)M\ (A.3) 

dtMi = -{l/2)a{din)-{dia)n + f3^ViMi + aKMi + {Vif3i)M^ 

= -{l/2)a{Vi'H) - {Via)H + P^ViMi + aKMi + iViPi)M\ (A.4) 

or by using Lie derivatives along an'^, 

Xant^n = 2aKn-2a-i-^l^di{^M^)-A{dia)M\ (A.5) 
£ani^Mi = -{l/2)a{diH)-{dia)H + aKMi. (A.6) 

Expression using jij and Kij In order to check the effects of the adjustments in (2.5) and (2.6) 
to constraint propagation, it is useful to re-express (A.l) and (A.2) using jij and Kij. By a straight- 
forward calculation, we obtain an expression as 

dtH = Hr{darnn) + Hr"di{darnn)+Hi^^^didj{darn^ (A.7) 
dtMi = Mii'^^idtJmn) + M2i^'^''dj{dt^mn) + Mu^^{dtKmn) + Mj'^^djidtKmn), (A.8) 

where 

jjmn 2^(3)"*" pP .p'^.^™*'y'^-? -|- p'^p'^ 

^Y^^^P^d^^mk^^Oj^^^) - 7™^™(5,7'^)(5,7fcp) - 2KK^- + 2K^jK^\ (A.9) 



^imn ._ _27™r"- (3/2)7^^(5^7"*")+ 7"^^' (6>j-y")+7"^"r, (A.IO) 

^jjmn _ _^ yn^mi^ (A.ll) 

/ff" := 2(^7™" - i^™"), (A.12) 

Mh™" := 7"'(5i^™j)-7"'(5ii^"i) + (V2)(5i7™")^'i + r"i^"i, (A.13) 
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M^Jmn ._ + (i/2)^™"Ki. + (l/2)K'""(5/, (A.14) 

Mai*"" := -5,"r"^-(l/2)(9i7""), (A.15) 
M4i^'"^" := 7™J(^^" -7™"^/, (A.16) 

where we expressed F"* = T^^''^ . 

A. 2 Constraint propagations for the adjusted ADM systems 

Generally, we here write the adjustment terms to (2.5) and (2.6) using (2.7) and (2.8) by the following 
combinations, [these are the same with (3.7) and (3.8)] 

adjustment term of dtjij : +Pijn + Q^^jMk + p''ij{Vk'H) + q''\j{VkMi), (A.17) 
adjustment term of dtKij : +Rij'H + S^ijMk + r''ij{VkH) + s'^^ijiVkMi), (A.18) 

where P, Q, R, S and p, q, r, s are multipliers (please do not confuse Rij with three Ricci curvature 

that we write as Rif)- We adjust them only using up to the first derivatives in order to make the 
discussion simple. 

By substituting the above adjustments into (A. 7) and (A. 8), we can write the adjusted constraint 
propagation equations as 

dtTC = (original terms) 

+Hr[Prnnn + Qtn^k + /™„(VfcW) + g'^'^n (VfcA^;)] 
+ifr"9.[PmnW + QL-^fc +/mn(Vfc7^) + (Z' mn(VfcA10] 

+Hr[Rmnn + St^Mk + r'=^n(VfcK) + s'='mn(VfeM,)], (A.19) 

dtM-i = (original terms) 

+Mir^[Pmnn + QinMh+p'^mni'^kH) + q'^^ni'^ kMl)] 
+M2i^'^^dj[Pmn'H + QtnMk +/mn(Vfe7^) + q^^mniVkMl)] 

+Msr''[Rmnn + Si^Mk + r'^mni^kn) + S^Ln(Vfc>lz)] 

+M^^'^''dj[Rmnn + Si^Mk + r'^mniVkn) + S*^'mn(VfeM,)]. (A.20) 

Here the "original terms" can be understood either as (A.l) and (A.2), or as (A.7) and (A. 8). There- 
fore, for example, we can see that adjustments to df-jij do not always keep the constraint propagation 
equations in the first order form, due to their contribution in the third adjusted term in (A.19). 

We note that these expressions of constraint propagation equations are equivalent when we include 
the cosmological constant and/or matter terms. 
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B Numerical demonstrations using the Ashtekar formulation 



This appendix is devoted to introduce our numerical comparisons between three levels of hyperbolicity 
using Ashtekar variables. Details are available in [77, 94]. 

B.l The Ashtekar formulation 

The key feature of Ashtekar's formulation of general relativity [13] is the introduction of a self-dual 
connection as one of the basic dynamical variables. Let us write the metric g^^ using the tetrad i?^ 
as g^y = E^^E^rjij ^. Define its inverse, by E'^ := E;;^g>^'^riij and we impose E^ = as the gauge 
condition. We define S0(3,C) connections ^yl^ := =F (i/2)e";,c'^)^'^) where ujj/ is a spin connection 
1-form (Ricci connection), ujj^^ := E^'^'V ^E^ . Ashtekar's plan is to use only the self-dual part of the 
connection ^A'^ and to use its spatial part '^Af as a dynamical variable. Hereafter, we simply denote 
'^Ajj^ as .4.^. 

The lapse function, N, and shift vector, N'^, both of which we treat as real- valued functions^, are 
expressed as Eq = {1/N, —N^/N). This allows us to think of Eq as a normal vector field to S spanned 
by the condition t = =const., which plays the same role as that of the ADM formulation. Ashtekar 
treated the set , AD as basic dynamical variables, where is an inverse of the densitized triad 
defined by Ej^ := e£^* , where e := det Ef is a density * . This pair forms the canonical set. In the case 
of pure gravitational spacetime, the Hilbert action takes the form 

S = j <^x\{dtAr)K + iimKHF-/'^c-N'F^^ + ^V,&^, (B.21) 

where N := e^'^N, F^^ := 25[^^^j - ie"-i,cA\Al is the curvature 2-form, ViE^ := diEi - ieab" A\Ei. 
The action (B.21) gives us the following evolution equations and constraints: 



The Ashtekar formulation [13]: Box B.l 

The dynamical variables are (£'*,^f). 

The evolution equations for a set of (£'^,.4") are 

dtK = -W,{e-\NMn) + 2'r>j{N^'&}) + iAUab''K, (B.22) 

dtAf = -ie'^'cNElF^j + N^F^. + ViA^, (B.23) 

where VjXl' := djXi' - ieab^A^Xl^ and F9-. ■- 2%^^"j - ie«6c A^-^j- 
Constraint equations: (Hamiltonian, momentum and Gauss constraints) 

Cf^ := m^'^'cKElF^j^^, (B.24) 

C^f ,= -Ft.Ei^O, (B.25) 

C^f := V.El^O, (B.26) 



®We use I,J = (0),---,(3) and a,b = (l),---,(3) arc 50(1,3), 50(3) indices respectively. Wc raise and lower 
fj,,v, - ■ ■ by g^" and g^i,i, (the Lorentzian metric); I, J, - ■ ■ by rj^' = diag(— 1, 1, 1, 1) and r/ij; ■ • ■ by 7*-' and jij (the 
three- metric); a, 6, ■ ■ ■ by (J"** and Sab- We also use volume forms eatc- ^atc^"''"^ ~ 3!. 

''N and A''' are the same with a and /J* in the previous text. We follow this conventional notation in this appendix. 

*For later convenience, = det.B* = {detEff = {l/^je"^" eijuKElE'^, where tijk ■= eabcE^E^E^ and €ijk := 

e~^e<3fc. When (i,j,k) = (1,2,3), we have eijk = e, eijk = 1, e'^'^ = e~^, and e*^* = 1. 
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Wc have to consider the reaUty conditions when we use this formahsm to describe the classical 
Lorentzian spacetime. Fortunately, the metric will remain on its real-valued constraint surface during 
time evolution automatically if we prepare initial data which satisfies the reality condition. More 
practically, we further require that triad is real-valued. But again this reality condition appears as 
a gauge restriction on ^q[91], which can be imposed at every time step. In our actual simulation, 
we prepare our initial data using the standard ADM approach, so that we have no difficulties in 
maintaining these reality conditions. 

B.2 Reformulate the Ashtekcir evolution equations 
B.2.1 Strongly and Symmetric Hyperbolic systems 

The authors' recent studies showed the following: 

(a) The original set of dynamical equations (B.22) and (B.23) [the original equations] already forms 
a weakly hyperbolic system [93]. So that we regard the mathematical structure of the original 
equations as one step advanced from the standard ADM. 

(b) Further, we can construct higher levels of hyperbolic systems by restricting the gauge condition 
and/or by adding constraint terms, C^^^,Cj^^ and Cq^^, to the original equations. 

— by requiring additional gauge conditions or adding constraints to the dynamical equations, 
we can obtain a strongly hyperbolic system [93], 

— by requiring additional gauge conditions and adding constraints to the dynamical equations, 
we can obtain a symmetric hyperbolic system [92, 93]. 

(c) Based on the above symmetric hyperbolic system, we can construct an Ashtckar version of the A- 
system [24] which is robust against perturbative errors for both constraints and reality conditions 
[76]. 

In order to obtain a symmetric hyperbolic system ^, we add constraint terms to the right-hand-side 
of (B.22) and (B.23). The adjusted dynamical equations, 

dtK = -iVj{e'\N&cM) + '^'I^jiN^E'I) + iA'oeab'K + ^^^^ (B.27) 

where P\b = N'Sab + iNeab^K^ 

dtAf = -ie''\NEiF^^ + N^Ff^ + ViA^ + K'iQ1Cf^ + K^Ri^''C^f, (B.28) 

where = e'^NE^, Ri^" = ie-^Ne''\E\Ei 

form a symmetric hyperbolicity if we further require ki = K2 = /^s = 1 and the gauge conditions, 

Al = AfN\ diN = 0. (B.29) 

We remark that the adjusted coefficients, P'''ab,Qi, Ri''"', for constructing the symmetric hyperbolic 
system are uniquely determined, and there are no other additional terms (say, no Cjf^,C^^ for dt&a^ 
no Cq^^ for dtAf) [93]. The gauge conditions, (B.29), are consequences of the consistency with (triad) 
reality conditions. 

We can also construct a strongly (or diagonalizable) hyperbolic system by restricting to a gauge 
0, ±N\/^ (where 7'' is the three-metric and we do not sum indices here) for the original 

^Iriondo et al [48] presented a symmetric hyperbolic expression in a different form. The differences between ours and 
theirs are discussed in [92, 93]. 
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system 


variables 


Eqs of motion 


remark 


I Ashtekar (weakly hyp.) 




(B.22), (B.23) (original) 


"original" eqs. 


II Ashtekar (strongly hyp.) 




(B.27), (B.28) (with k=1) 


(B.30) required 


Ill Ashtekar (symmetric hyp.) 




(B.27), (B.28) (with k = 1) 


(B.29) required 


adj Ashtekar (adjusted) 




(B.27), (B.28) (with k ^ 1) 




A Ashtekar-A-system 




(B.37) 


controls CH,CMi,CGa 



Table 5: List of systems that we compare in this appendix. 



equations (B.22), (B.23). Or we can also construct from the adjusted equations, (B.27) and (B.28), 
together with the gauge condition 



As for the strongly hyperbolic system, we hereafter take the latter expression. 
In Table 5, we summarized the equations to be used for our comparisons. 



(B.30) 



B.2.2 Ashtekar-A-system 

In §2.3.1, we introduced an idea to construct a robust evolution system against a perturbative error, 
named "A-systcm'" [24]. Based on the above symmetric hyperbolic equations, we constructed the 
Ashtekar version of the "A-system" [76]. Here we present our system which evolves the spacetime to 
the constraint surface, Ch ~ Cmi ~ Cca ~ as the attractor. In [76], we also presented a system 
which controls the perturbative violation of the reality condition. 

We introduce new variables (A, Aj, Aq), as they obey the dissipative evolution equations 

dt\ = aiCH-PiX, (B.31) 
dtXi = a2CMi-p2Xu (B.32) 
dtK = asCca-Z^sAa, (B.33) 

where ai ^ (allowed to be complex numbers) and Pi > (real numbers) are constants. 



If we take u 



{DL) 



{El,Af,X,Xi,Xa) as a set of dynamical variables, then the principal part of 



(B.31)-(B.33) can be written as 



dtX 
dtXi 
dtXa 



J)cd i 



taie--E^^E'a{diAj), (B.34) 
a2[-e6iEl+e6iEi]{diA'j), (B.35) 
a3diEl (B.36) 

The characteristic matrix of the system u^^'' does not form a Hermitian matrix. However, if we 

modify the right-hand-side of the evolution equation of {El^,Af), then the set becomes a symmet- 
ric hyperbolic system. This is done by adding a^Y^diXa) to the equation of dtE^^, and by adding 
iaie^/EI&^idiX) + a2(-ey™4" + edf &''){diXm) to the equation of dtA^. The final principal part, 
then, is written as 

-iaieb'"^E^&^ 
a2e{SrEi - 6\E^) 
Ka^daH^m y 



a, 



A 
Ai 

Va„/ 



^ 






















Ab \ 



di 



A 

Am 

Va, / 

(B.37) 
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Clearly, the solution {E^, , X, Xi, Xa) = {E^, Af ,0,0,0) represents the original solution of the 
Ashtekar system. If the As decay to zero after the evolution, then the solution also describes the 
original solution of the Ashtekar system in that stage. Since the dynamical system of u^^\ (B.37), 
constitutes a symmetric hyperbolic form, the solutions to the A-system are unique. Therefore, the 
dynamical system, (B.37), is useful for stabilizing numerical simulations from the point that it recovers 
the constraint surface automatically. 

B.2.3 Adjusted Ashtekar system 

We also try to compare a set of evolution system, which we proposed as "adjusted-system". The 
fundamental equations that we will demonstrate are the same with (B.27) and (B.28), but here the 
real-valued constant multipliers ks are not necessary equals to unity. We set k = ki = K2 = for 
simplicity. Apparently the set of (B.27) and (B.28) becomes the original weakly hyperbolic system if 
K = 0, becomes the symmetric hyperbolic system if k = 1 and N = const.. The set remains strongly 
hyperbolic systems for other choices of k except k = 1/2 which only forms weakly hyperbolic system. 

B.3 Comparing numerical performance 
B.3.1 Model and Numerical method 

The model we present here is gravitational wave propagation in a planar spacetime under periodic 
boundary condition. Wc performed a full numerical simulation using Ashtckar's variables. We prepare 
two +- mode strong pulse waves initially by solving the ADM Hamiltonian constraint equation, using 
York-O'Murchadha's conformal approach. Then we transform the initial Cauchy data (3-metric and 
extrinsic curvature) into the connection variables, (-£^,.4"), and evolve them using the dynamical 
equations. For the presentation in this article, we apply the geodesic slicing condition (ADM lapse 
A'' = 1 or densitized lapse N = 1, with zero shift and zero triad lapse). We have used both the 
Brailovskaya integration scheme, which is a second order predictor-corrector method, and the so- 
called iterative Crank-Nicholson integration scheme for numerical time evolution. The details of the 
numerical method arc described in [77]. 

More specifically, we set our initial guess 3-metric as 



in the periodically bounded region x = [—5, -|-5]. Here K and L are constants and we set K = 0.3 and 
L = 2.5 for the plots. 

In order to show the expected "stabilization behavior" clearly, we artificially add an error in the 
middle of the time evolution. We added an artificial inconsistent rescaling once at time t = 6 for the 
Ay component as Ay — > Ay{l + error). 

B.3. 2 Differences between three levels of hyperbolicity 

We have performed comparisons of stability and/or accuracy between weakly and strongly hyperbolic 
systems, and between weakly and symmetric hyperbolic systems[77]. (We can not compare strongly 
and symmetric hyperbolic systems directly, because these two requires different gauge conditions.) 

We omit figures in this report, but one can see a part of results in Fig.9 and Fig. 10. We may 
conclude that higher level hyperbolic system gives us slightly accurate evolution. However, if we 
evaluate the magnitude of L2 norms, then we also conclude that there is no measurable differences 
between strongly and symmetric hyperbolicities. This last fact will be supported more affirmatively 
in the next experiments. 






l + i^(e-(^-'^)' +e-(^+-^)') 




(B.38) 



sym. 
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Figure 9: Demonstration of the Ashtekar-A-system for the cases of plane wave propagation under the periodic 
boundary. We plot the L2 norm of the Hamiltonian constraint equation, Ch- Fig. (a) shows how the system 
goes bad depending on the amplitude of artificial error at t = 6. All the lines are of the evolution of Ashtekar's 
original equation (no A-system). Fig. (b) shows the effect of A-system. All the lines are 20% error amplitude, 
but shows the difference of evolution equations. The solid line is for Ashtekar's original equation (the same as in 
Fig. (a)), the dotted line is for the strongly hyperbolic Ashtekar's equation. Other lines are of A-systems, which 
produces better performance than that of the strongly hyperbolic system. (Reprinted from [94], ©lOP 2001) 



B.3.3 Demonstrating "A-system" 

Next, we show a result of the "A-system" [94]. Fig. 9 (a) shows how the violation of the Hamiltonian 
constraint equation, Ch, become worse depending on the term error. The oscillation of the L2 norm 
Ch in the figure due to the pulse waves collide periodically in the numerical region. We, then, fix the 
error term as a 20% spike, and try to evolve the same data in different equations of motion, i.e., the 
original Ashtekar's equation [solid line in Fig. 9 (b)], strongly hyperbolic version of Ashtekar's equation 
(dotted line) and the above A-system equation (other lines) with different (3s but the same a. As we 
expected, all the A-system cases result in reducing the Hamiltonian constraint errors. 




Figure 10: Demonstration of the adjusted-Ashtekar systems. The experiments as in Fig. 9(b). Fig. (a) and (b) 
are L2 norm of the Hamiltonian constraint equation, Ch, and momentum constraint equation, Cmx, respectively. 
The solid line is the case of k = 0, that is the case of "no adjusted" original Ashtekar equation (weakly hyperbolic 
system). The dotted line is for k = 1, equivalent to the symmetric hyperbolic system. We see other line (k = 2.0) 
shows better performance than the symmetric hyperbolic case. (Reprinted from [94], ©lOP 2001) 
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B.3.4 Demonstrating "adjusted system" 



We here examine how the adjusted multipHers contribute to the system's stabiHty [94]. In Fig. 10, we 
show the results of this experiment. We plot the violation of the constraint equations both Ch and 
Cmx- An artificial error term was added in the same way as above. The solid line is the case of «; = 0, 
that is the case of "no adjusted" original Ashtekar equation (weakly hyperbolic system). The dotted 
line is for n = I, equivalent to the symmetric hyperbolic system. We see other line (k = 2.0) shows 
better performance than the symmetric hyperbolic case. 
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